[R] integrate function - error -integration not occurring with last few rows

Berend Hasselman bhh at xs4all.nl
Fri Apr 6 07:47:41 CEST 2012


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



More information about the R-help mailing list