Creation of Animal category in p.d solved all problems. Plots fine now. The smallest hurdles are often the hardest to get over.
Gavin Simpson wrote: > > On Mon, 2009-05-18 at 11:48 -0700, William Paterson wrote: >> Hi, >> >> I am using GAMMs to show a relationship of temperature differential over >> time with a model that looks like this:- >> >> gamm(Diff~s(DaysPT)+AirToC,method="REML") >> >> where DaysPT is time in days since injury and Diff is repeat measures of >> temperature differentials with regards to injury sites compared to >> non-injured sites in individuals over the course of 0-24 days. I use the >> following code to plot this model on the response scale with 95% CIs >> which >> works fine:- >> >> g.m<-gamm(Diff~s(DaysPT)+AirToC,method="REML") >> p.d<-data.frame(DaysPT=seq(min(DaysPT),max(DaysPT))) >> p.d$AirToC<-(6.7) >> b<-predict.gam(g.m$gam,p.d,se=TRUE) >> range<-c(min(b$fit-2*b$se.fit),max(b$fit+2*b$se.fit)) >> plot(p.d$DaysPT,b$fit,ylim=c(-4,12),xlab="Days post-tagging",ylab="dTmax >> (ºC)",type="l",lab=c(24,4,12),las=1,cex.lab=1.5, cex.axis=1,lwd=2) >> lines(p.d$DaysPT,b$fit+b$se.fit*1.96,lty=2,lwd=1.5) >> lines(p.d$DaysPT,b$fit-b$se.fit*1.96,lty=2,lwd=1.5) >> points(DaysPT,Diff) >> >> >> However, when I add a correlation structure and/or a variance structure >> so >> that the model may look like:- >> >> >> gamm(Diff~s(DaysPT3)+AirToC,correlation=corCAR1(form=~DaysPT| >> Animal),weights=varPower(form=~DaysPT),method="REML") >> >> >> I get this message at the point of inputting the line >> "b<-predict.gam(g.m$gam,p.d,se=TRUE)" > > Note that p.d doesn't contain Animal. Not sure this is the problem, but > I would have thought you'd need to supply new values of Animal for the > data you wish to predict for in order to get the CAR(1) errors correct. > Is it possible that the model is finding another Animal variable in the > global environment? > > I have predicted from several thousand GAMMs containing correlation > structures similar to the way you do above so this does work in general. > If the above change to p.d doesn't work, you'll probably need to speak > to Simon Wood to take this further. > > Is mgcv up-to-date? I am using 1.5-5 that was released in the last week > or so. > > For example, this dummy example runs without error for me and is similar > to your model > > y1 <- arima.sim(list(order = c(1,0,0), ar = 0.5), n = 200, sd = 1) > y2 <- arima.sim(list(order = c(1,0,0), ar = 0.8), n = 200, sd = 3) > x1 <- rnorm(200) > x2 <- rnorm(200) > ind <- rep(1:2, each = 200) > d <- data.frame(Y = c(y1,y2), X = c(x1,x2), ind = ind, time = rep(1:200, > times = 2)) > require(mgcv) > mod <- gamm(Y ~ s(X), data = d, corr = corCAR1(form = ~ time | ind), > weights = varPower(form = ~ time)) > p.d <- data.frame(X = rep(seq(min(d$X), max(d$X), len = 20), 2), > ind = rep(1:2, each = 20), > time = rep(1:20, times = 2)) > pred <- predict(mod$gam, newdata = p.d, se = TRUE) > > Does this work for you? If not, the above would be a reproducible > example (as asked for in the posting guide) and might help Simon track > down the problem if you are running an up-to-date mgcv. > > HTH > > G > >> >> >> Error in model.frame(formula, rownames, variables, varnames, extras, >> extranames, : >> variable lengths differ (found for 'DaysPT') >> In addition: Warning messages: >> 1: not all required variables have been supplied in newdata! >> in: predict.gam(g.m$gam, p.d, se = TRUE) >> 2: 'newdata' had 25 rows but variable(s) found have 248 rows >> >> >> Is it possible to predict a more complicated model like this on the >> response >> scale? How can I incorporate a correlation structure and variance >> structure >> in a dataframe when using the predict function for GAMMs? >> >> Any help would be greatly appreciated. >> >> William Paterson >> >> >> >> > -- > %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~% > Dr. Gavin Simpson [t] +44 (0)20 7679 0522 > ECRC, UCL Geography, [f] +44 (0)20 7679 0565 > Pearson Building, [e] gavin.simpsonATNOSPAMucl.ac.uk > Gower Street, London [w] http://www.ucl.ac.uk/~ucfagls/ > UK. WC1E 6BT. [w] http://www.freshwaters.org.uk > %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~% > > ______________________________________________ > 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. > > -- View this message in context: http://www.nabble.com/Predicting-complicated-GAMMs-on-response-scale-tp23603248p23639184.html Sent from the R help mailing list archive at Nabble.com. ______________________________________________ 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.
