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.
