Apologies for the lengthy code. I tried a simple (and shorter) piece of code (pasted below) and it still gives me the same error for last few rows. Is this a bug or am I doing something totally wrong? Could anyone please provide some help/pointers ?
PS. beta0 was fixed to 0.001 in the previous code. Also if I manually estimated the integral isnt 0. If I scramble the row order, it is again only the last few rows that dont integrate. Thanks data1<-expand.grid(1:5,0:6,10) names(data1)<-c("ID","TIME", "DOSE") data1<-data1[order(data1$DOSE,data1$ID,data1$TIME),] ################ ed<-data1[!duplicated(data1$ID) , c(1,3)] ed$base=1 ed$drop=1 set.seed(5234123) k<-0 for (i in 1:length(ed$ID)) { k<-k+1 ed$base[k]<-100*exp(rnorm(1,0,0.2)) ed$drop[k]<-0.20*exp(rnorm(1,0,0.5)) } comb1<-merge(data1[, c("ID","TIME")], ed) comb1$disprog<-comb1$base*exp(-comb1$drop*comb1$TIME) comb1$integral=0 hz.func1<-function(t,bshz,beta1, change) { ifelse(t==0,bshz, bshz*exp(beta1*change)) } q<-0 for (m in i:length(comb1$ID)) { q<-q+1 comb1$integral[q]<-integrate(hz.func1, lower=0, upper=comb1$TIME[q], bshz=0.001,beta1=0.035, change=comb1$disprog[q])$value } comb1 On Fri, Apr 6, 2012 at 1:47 AM, Berend Hasselman <b...@xs4all.nl> wrote: > > 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 > > -- 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.