Hi Tim,

How does this sound for the email to the list:




Dear list users,





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]]


