From: Simon Wood <s.wood_at_bath.ac.uk>

Date: Mon, 27 Jun 2011 09:49:22 +0100

R-help_at_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. Received on Mon 27 Jun 2011 - 08:51:44 GMT

Date: Mon, 27 Jun 2011 09:49:22 +0100

> In your dummy variable trick, I see how subjectwise predictions are

*> obtained. But to get the Group mean predictions, isn't the s term
**> with the 0 dummy by variable the same as just omitting the s term
**> with Subject altogether? I would think so but want to check.
*

- Yes it's the same as omitting the term altogether for prediction (of
means), but that's not the same as having omitted it for fitting, of
course.

> In the spirit of handling subjectwise variabilities, I could imagine

*> using the term
**>
**> s(Subject, bs="re", by=X)
**>
**> and trying to emulate lme-like behavior and fitting a random effect
**> slope to each subject. First question is, is that indeed what it
**> does? If so, would this term include an intercept? I would assume
**> it doesn't but again want to check.
*

- s(Subject,X,bs="re") is a slightly better way of doing this, as it
generalises a bit more easily. To include random intercepts as well use

s(Subject,bs="re") + s(Subject,X,bs="re") ...exactly how these terms work is explained in the details section of ?smooth.construct.re.smooth.spec.

best,

Simon

> Thanks again, John

*>
**>
**>
**> -----Original Message----- From: Simon Wood
**> [mailto:s.wood_at_bath.ac.uk] Sent: Friday, 24 June, 2011 11:43 AM To:
**> Szumiloski, John Cc: r-help_at_r-project.org Subject: Re: mgcv:gamm:
**> predict to reflect random s() effects?
**>
**> Given that you don't have huge numbers of subjects you could fit the
**> model with `gam' rather than `gamm', using
**>
**> out.gamm<- gam( Y ~ Group + s(X, by=Group) + s(Subject,bs="re"),
**> method="REML")
**>
**> Then your predictions will differ by subject (see e.g.
**> ?random.effects for a bit more information on simple random effects
**> in mgcv:gam).
**>
**> A further trick allows you to choose whether to predict with the
**> subject effects at their predicted values, or zero.
**>
**> Let dum be a vector of 1's...
**>
**> out.gamm<- gam( Y ~ Group + s(X, by=Group) +
**> s(Subject,bs="re",by=dum), method="REML")
**>
**> Predicting with dum set to 1 gives the predictions that you want.
**> Setting dum to 0 gives predictions with the prediction Subject
**> effects set to zero.
**>
**> The reason that trying to predict with the gamm lme object is tricky
**> relates to how gamm works. It takes the GAMM specification, and then
**> sets up a corresponding `working mixed model' which is estimated
**> using lme. The working mixed model uses working variable names set
**> within gamm. If you try to predict using the working model lme
**> object then predict.lme looks for these internal working variable
**> names, not the variable names that you supplied....
**>
**> Basically gamm treats all random effects as 'part of the noise' in
**> the model specification, and adjusts the variance estimates for the
**> smooths and fixed effects to reflect this. It isn't set up to
**> predict easily at different random effect grouping levels, in the way
**> that lme is.
**>
**> best, Simon
**>
**> On 24/06/11 15:59, Szumiloski, John wrote:
**>> Dear useRs, I am using the gamm function in the mgcv package to
**>> model a smooth relationship between a covariate and my dependent
**>> variable, while allowing for quantification of the subjectwise
**>> variability in the smooths. What I would like to do is to make
**>> subjectwise predictions for plotting purposes which account for
**>> the random smooth components of the fit. An example. (sessionInfo()
**>> is at bottom of message) My model is analogous to> out.gamm<-
**>> gamm( Y ~ Group + s(X, by=Group), random = list(Subject=~1) ) Y and
**>> X are numeric, Group is an unordered factor with 5 levels, and
**>> Subject is an unordered factor with ~70 levels Now the output from
**>> gamm is a list with an lme component and a gam component. If I make
**>> a data frame "newdat" like this:
**>>> newdat
**>> X Group Subject 5 g1 s1 5 g1 s2 5 g1 s3 6 g1 s1 6 g1 s2 6 g1 s3 I
**>> can get the fixed effects prediction of the smooth by>
**>> predict(out.gamm$gam, newdata=newdat) Which gives 1 1.1 1.2 2 2.1
**>> 2.2 3.573210 3.573210 3.573210 3.553694 3.553694 3.553694 But I
**>> note that the predictions are identical across different values of
**>> Subject. So this accounts for only the fixed effects part of the
**>> model, and not any random smooth effects. If I try to extract
**>> predictions from the lme component:
**>>> predict(out.gamm$lme, newdata=newdat) I get the following error
**>> message: Error in predict.lme(out.gamm$lme, newdata = newdat) :
**>> Cannot evaluate groups for desired levels on "newdata" So, is
**>> there a way to get subjectwise predictions which include the
**>> random effect contributions of the smooths? Thanks, John ---------
**>> ### session info follows
**>>> sessionInfo()
**>> R version 2.13.0 Patched (2011-06-20 r56188) Platform:
**>> i386-pc-mingw32/i386 (32-bit) locale: [1]
**>> LC_COLLATE=English_United States.1252 LC_CTYPE=English_United
**>> States.1252 [3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C
**>> [5] LC_TIME=English_United States.1252 attached base packages: [1]
**>> grDevices datasets splines grid graphics utils stats methods base
**>> other attached packages: [1] mgcv_1.7-6 gmodels_2.15.1 car_2.0-10
**>> nnet_7.3-1 MASS_7.3-13 nlme_3.1-101 [7] rms_3.3-1 Hmisc_3.8-3
**>> survival_2.36-9 lattice_0.19-26 loaded via a namespace (and not
**>> attached): [1] cluster_1.14.0 gdata_2.8.2 gtools_2.6.2
**>> Matrix_0.999375-50 tools_2.13.0 John Szumiloski, Ph.D.
**> [snip]
**>
**>
**> -- Simon Wood, Mathematical Science, University of Bath BA2 7AY UK
**> +44 (0)1225 386603 http://people.bath.ac.uk/sw283
**> Notice: This e-mail message, together with any attach...{{dropped:18}}
*

R-help_at_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. Received on Mon 27 Jun 2011 - 08:51:44 GMT

This quarter's messages: by month, or sorted: [ by date ] [ by thread ] [ by subject ] [ by author ]

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.2.0, at Mon 27 Jun 2011 - 10:00:17 GMT.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help.
Please read the posting
guide before posting to the list.
*