[R-sig-ME] Numerically integrating a GAMM from expected size-specific growth to get size-at-age growth curves

[Nicholas Lewin-Koh]

Hi Rebecca,
I am not familiar with these particular papers but we have used gamms
to model tumor growth in tumor xenograft models, where our collaborators
are interested
in the AUC of the growth curve. If all you are interested in
are the fixed effects, average growth in the population, it is not too
hard. If you look at the help page for predict.gam Simon has a nice
example
doing something similar to what you are doing, including getting a
confidence interval
for the curve.

I don't understand in your code why you need to add eps? In the gam part
of
the model fit the variation from the random effects is incorporated
in the estimates, Wouldn't something like;

hl1<-gamm(ahlgr~s(hl),family= gaussian, data=ahl,random=list(id=~1))
pred1<-predict(hl1$gam,se=T,type="response")
m<-which(pred1==median(pred1))
trapInt(pred1[1:m], delta=.01) ## using trapezoidal rule for
integration, this is for the mean, up to the median growth assuming your
data is sorted
x.mesh <- seq(0,.5,length=100)
newd <- data.frame(hl = x.mesh)
library(mvtnorm)
br <- rmvnorm(1000,coef(hl1$gam),hl1$gam$Vp)
pp<-2^(Xp%*%t(br))
simInt <- apply(pp,2,function(y)tapply(y,trapInt,delta=.01))
lowup <- round(t(apply(simInt,1,quantile,probs =
c(.025,1-.025),names=FALSE))) 

The covariance is estimated on the scale of the data, and is not 
correct for the integrated response, hence the need to simulate, or do 
the nasty mathe to get the transformed covariance using the delta
method.

Hope this helps.

Nicholas


> Date: Sun, 4 Sep 2011 11:41:05 +1000
> From: Rebecca Laver <r.laver at student.unimelb.edu.au>
> To: r-sig-mixed-models at r-project.org
> Subject: [R-sig-ME] Numerically integrating a GAMM from expected
> 	size-specific growth to get size-at-age growth curves
> Message-ID:
> 	<CABwRE9LUV9z4okfudK91LOnTxHtGZa=fY9nKa1VGa81pNvJ_OA at mail.gmail.com>
> Content-Type: text/plain
> 
> Dear list users,
> 
> 
> 
> 
> 
> This has been sent out earlier but may have been addressed incorrectly,
> so
> here's a second try. I am an MSc  student using GAMMs for analysing
> somatic
> growth rate patterns in Komodo dragons (using a reasonably large data set
> that my supervisor and his colleagues have collected over the last 9
> years)
> in an  attempt to understand how spatial  and individual based covariates
> could influence key life-history parameters in these lizards. Similar
> statistical approaches have been used elsewhere (eg. sea turtles studies
> analysed by Milani Chaloupka)  and have provided an important way of
> modelling growth in animals that have more complex growth trajectories
> than
> predicted by classic growth models  (eg. von Bertalanffy and Richard’s
> models). We would really like to be able to numerically integrate these
> models to predict growth trajectories with respect to time based
> estimates
> (i.e. how many years does it take to reach maturity or maximal body size;
> such as those in Fig. 8 b and c in Balazs, G. H. and Chaloupka, M., 2004
> ~
> http://www.seaturtle.org/PDF/Balazs_2004_MarBiol.pdf). It is likely to
> help
> us understand, what we expect is a huge degree of individual and spatial
> variation in growth trajectories for these large predatory lizards.
> 
> 
> 
> I realise that a very similar topic  came up on the list last year. We
> have
> followed up on these comments with trying to use them to work this issue
> out,  and we even contacted the subscriber who requested help, and she
> also
> reported that after a couple of weeks of trying that it did not work out.
> So
> if I can please again request help with this issue that would be great.
> 
> 
> 
> Sincerely,
> 
> Rebecca
> 
> 
> P.s. here is the code we are using and we have tried for the GAMM.
> 
> #Creating GAMM model and plotting on original scale#
> hl1<-gamm(ahlgr~s(hl),family=
> gaussian, data=ahl,random=list(id=~1))
> pred1<-predict(hl1$gam,se=T,type="response")
> l1<-order(ahl$hl)
> plot(ahl$hl,ahlgr,cex=0.5,col="dark grey",xlab="Head Length (cm)",
> ylab="Growth Rate (cm HL/yr)")
> lines(ahl$hl[l1],pred1$fit[l1],lwd=2)
> lines(ahl$hl[l1],pred1$fit[l1]+2*pred1$se.fit[l1],lty="dashed")
> lines(ahl$hl[l1],pred1$fit[l1]-2*pred1$se.fit[l1],lty="dashed")
> 
> #We tried the following after this was suggested by Simon Wood in
> response
> to another subscriber's previous request for assistance: "It's best to do
> the integration and differentiation numerically. See the example at the
> end
> of ?predict.gam, for a differentiation example. To integrate I usually
> use
> predict.gam to evaluate the smooth on a fine mesh, and then just sum the
> values and multiply by the mesh spacing (this is the midpoint rule)"#
> 
> x.mesh <- seq(0,1,length=848)
> newd <- data.frame(hl = x.mesh)
> X0 <- predict(hl1$gam,newd,type="lpmatrix")
> eps <- 1e-7
> x.mesh <- x.mesh + eps
> newd <- data.frame(hl = x.mesh)
> X1 <- predict(hl1$gam,newd,type="lpmatrix")
> Xp <- (X1-X0)/eps
> colnames(Xp)
> for (i in 1) {
>   Xi <- Xp*0
>    Xi[,(i-1)*9+1:9+1] <- Xp[,(i-1)*9+1:9+1]
>     df <- Xi%*%coef(hl1$gam)
>      df.sd <- rowSums(Xi%*%hl1$gam$Vp*Xi)^.5
>       plot(x.mesh,df,type="l",ylim=range(c(df+2*df.sd,df-2*df.sd)))
>        lines(x.mesh,df+2*df.sd,lty=2);lines(x.mesh,df-2*df.sd,lty=2)
>         }
> 
> #However, as I have had some difficulty understanding the process of
> evaluating smooths on a fine mesh I am uncertain as to whether I have
> gone
> through further steps in the example than I need to, and am unsure what
> to
> do next. #
> 
> 	[[alternative HTML version deleted]]
> 
> 
> 
> ------------------------------
> 
> _______________________________________________
> R-sig-mixed-models mailing list
> R-sig-mixed-models at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> 
> 
> End of R-sig-mixed-models Digest, Vol 57, Issue 2
> *************************************************
>