On 06-04-2012, at 00:55, Navin Goyal wrote:
> 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
> }
Where is beta0 in the line with ed$bshz<-beta0 ?
In the last for loop for (m in i:length(comb3$ID)) could it be that i should
be 1?
When the i is changed to 1 then integrate results are <> 0, which might be what
you expect?
Berend
______________________________________________
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.
