Hi,
I am using the integrate function in some simulations in R (tried ver 2.12
and 2.15). The problem I have is that the last few rows do not integrate
correctly. I have pasted the code I used.
The column named "integral" shows the output from the integrate function.
The last few rows have no integration results. I tried increasing the
doses, number of subjects, etc.... this error occurs with the last few rows
only
I am not sure why this is happening. Could someone please help me with this
issue ??
Thank you for your time
dose<-c(0)
time<-(0:6)
id<-1:25
data1<-expand.grid(id,time,dose)
names(data1)<-c("ID","TIME", "DOSE")
data1<-data1[order(data1$DOSE,data1$ID,data1$TIME),]
################
basescore=95
basescore_sd=0.12
fall=0.15
fall_sd=0.5
slope=5
dose_slope1=0.045
dose_slope2=0.045
dose_slope3=0.002
rise_sd=0.5
ed<-data1[!duplicated(data1$ID) , c(1,3)]
ed$base=1
ed$drop=1
ed$bshz<-1
ed$up<-1
ed
set.seed(5234123)
k<-0
for (i in 1:length(ed$ID))
{
k<-k+1
ed$base[k]<-basescore*exp(rnorm(1,0,basescore_sd))
ed$drop[k]<-fall*exp(rnorm(1,0,fall_sd))
ed$up[k]<-slope*exp(rnorm(1,0,rise_sd))
ed$bshz<-beta0
}
comb1<-merge(data1[, c("ID","TIME")], ed)
comb1$disprog<-1
comb1$beta1<-0.035
comb1$beta21<-0.02
comb1$beta22<-0.45
comb1$beta23<-0085
comb1$beta31<-0.7
comb1$beta32<-0.05
comb1$exphz<-1
comb2<-comb1
p<-0
for(l in 1:length(comb2$ID))
{
p<-p+1
comb2$disprog[p]<-comb2$base[p]*exp(-comb2$drop[p]*comb2$TIME[p]) +
comb2$up[p]*comb2$TIME[p]
comb2$frac[p]<-ifelse ( comb2$DOSE[p]==3,
comb2$beta31[p]*comb2$TIME[p]^comb2$beta32[p],
exp(-comb2$beta21[p]*comb2$DOSE[p])*comb2$TIME[p]^comb2$beta22[p] )
}
hz.func1<-function(t,bshz,beta1, change,other)
{
ifelse(t==0,bshz, bshz*exp(beta1*change+other))
}
comb3<-comb2
comb3$integral=0
q<-0
for (m in i:length(comb3$ID))
{
q<-q+1
comb3$integral[q]<-integrate(hz.func1, lower=0, upper=comb3$TIME[q],
bshz=comb3$bshz[q],beta1=comb3$beta1[q],
change=comb3$disprog[q], other=comb3$frac[q])$value
}
#comb3[comb3$TIME==3, ] #
#tail(comb3)
--
Navin Goyal
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
