[R] Predicting complicated GAMMs on response scale
wdp1 at st-andrews.ac.uk
Wed May 20 18:40:54 CEST 2009
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:
>> I am using GAMMs to show a relationship of temperature differential over
>> time with a model that looks like this:-
>> 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
>> works fine:-
>> 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)
>> However, when I add a correlation structure and/or a variance structure
>> that the model may look like:-
>> I get this message at the point of inputting the line
> 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))
> 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.
>> 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
>> scale? How can I incorporate a correlation structure and variance
>> 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 at r-project.org mailing list
> PLEASE do read the posting guide
> 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.
More information about the R-help