Melinda Power
2011-Jul-19 21:16 UTC
[R] Incorrect degrees of freedom for splines using GAMM4?
Hello, I'm running mixed models in GAMM4 with 2 (non-nested) random intercepts and I want to include a spline term for one of my exposure variables. However, when I include a spline term, I always get reported degrees of freedom of less than 1, even when I know that my spline is using more than 1 degree of freedom. For example, here is the code for my model:> global.gamm4<-gamm4(zcog~s(adjpatx, fx=TRUE, k=5)+int234+cogagec+cogagesq+ + + oldfran +newus +alc2 +alc3 +alc4 +alcmiss +smk2 +smk3 + +mdinc10c +mdinc10sq+ pwhtc +pwhtsq +edu2+ edu3 +husbgs +husbcol+ husbmiss + +currpmh +pastpmh +neverpmh, random= ~(1|id) +(1|cogtest), data=global) Using > summary(global.gamm4$mer), I get the following output for my spline term, indicating that I use the expected 4 degrees of freedom. Xs(adjpatx)Fx1 0.1018943 0.1073225 0.949 Xs(adjpatx)Fx2 -0.0708114 0.1123845 -0.630 Xs(adjpatx)Fx3 0.7459511 0.6836413 1.091 Xs(adjpatx)Fx4 -0.2062321 0.0923569 -2.233 However, when I use > summary(global.gamm4$gam). I get an estimate of degrees of freedom that is not 4: Approximate significance of smooth terms: edf Ref.df F p-value s(adjpatx) 0.7588 0.7588 1.346 0.234 This degree of freedom = 0.76 also shows up on my plot. Ultimately, I would like to use a cubic regression penalized spline, allowing R to choose the degrees of freedom for me using GCV. However, when I use the correct code for this or variants of it using mgcv, I also get degrees of freedom less than 1. For example, in the following code provides a degree of freedom of less than 1 as well:> global.gamm4<-gamm4(zcog~s(adjpatx, fx=FALSE)+int234+cogagec+cogagesq ++ + oldfran +newus +alc2 +alc3 +alc4 +alcmiss +smk2 +smk3 + +mdinc10c +mdinc10sq+ pwhtc +pwhtsq +edu2+ edu3 +husbgs +husbcol+ husbmiss + +currpmh +pastpmh +neverpmh, random= ~(1|id) +(1|cogtest), data=global) Output indicating that this spline should probably look linear:> summary(global.gamm4$mer)Random effects: Groups Name Variance Std.Dev. id (Intercept) 0.1823454 0.427019 cogtest (Intercept) 0.0025498 0.050496 Xr.1 s(adjtibx) 0.0000000 0.000000 Residual 0.7782969 0.882211 Xs(adjtibx)Fx1 -0.0387360 0.0215596 -1.797 Output getting a df for this spline of 0.20.> summary(global.gamm4$gam)Approximate significance of smooth terms: edf Ref.df F p-value s(adjtibx) 0.2009 0.2009 16.07 NA The plot looks linear, but reports a df =0.20. So...to summarize my questions: 1. Are the splines produced by s(exp, fx=FALSE) or s(exp, fx=TRUE, k=k) correct even though the reported degrees of freedom appears to be wrong? 2. Can I believe my plot? 3. How can I get the true df used when I use s(exp, fx=FALSE)? Thanks for any and all help you can provide! Melinda [[alternative HTML version deleted]]
Simon Wood
2011-Jul-20 09:17 UTC
[R] Incorrect degrees of freedom for splines using GAMM4?
Melinda, Any chance you could send me the data that goes with this, offline, (if you can, I'll only use the data for investigating this issue, of course)? I can reproduce something going wrong with the edf computation for fixed smooths in gamm4, when there are random effects, but I can't get anywhere near as extreme a discrepancy as you are seeing. Also, I can't get anything to go wrong in mgcv, (actually you didn't paste the mgcv code, but rather another gamm4 call? could you send the mgcv code too?). Also, could you send me gamm4, mgcv and R version numbers, please? best, Simon On 19/07/11 22:16, Melinda Power wrote:> Hello, > > I'm running mixed models in GAMM4 with 2 (non-nested) random intercepts and > I want to include a spline term for one of my exposure variables. However, > when I include a spline term, I always get reported degrees of freedom of > less than 1, even when I know that my spline is using more than 1 degree of > freedom. For example, here is the code for my model: > >> global.gamm4<-gamm4(zcog~s(adjpatx, fx=TRUE, k=5)+int234+cogagec+cogagesq > + > + + oldfran +newus +alc2 +alc3 +alc4 +alcmiss +smk2 +smk3 > + +mdinc10c +mdinc10sq+ pwhtc +pwhtsq +edu2+ edu3 +husbgs > +husbcol+ husbmiss > + +currpmh +pastpmh +neverpmh, random= ~(1|id) > +(1|cogtest), data=global) > > Using> summary(global.gamm4$mer), I get the following output for my spline > term, indicating that I use the expected 4 degrees of freedom. > > Xs(adjpatx)Fx1 0.1018943 0.1073225 0.949 > Xs(adjpatx)Fx2 -0.0708114 0.1123845 -0.630 > Xs(adjpatx)Fx3 0.7459511 0.6836413 1.091 > Xs(adjpatx)Fx4 -0.2062321 0.0923569 -2.233 > > However, when I use> summary(global.gamm4$gam). I get an estimate of > degrees of freedom that is not 4: > > Approximate significance of smooth terms: > edf Ref.df F p-value > s(adjpatx) 0.7588 0.7588 1.346 0.234 > > This degree of freedom = 0.76 also shows up on my plot. > > Ultimately, I would like to use a cubic regression penalized spline, > allowing R to choose the degrees of freedom for me using GCV. However, when > I use the correct code for this or variants of it using mgcv, I also get > degrees of freedom less than 1. For example, in the following code provides > a degree of freedom of less than 1 as well: > > >> global.gamm4<-gamm4(zcog~s(adjpatx, fx=FALSE)+int234+cogagec+cogagesq + > + + oldfran +newus +alc2 +alc3 +alc4 +alcmiss +smk2 +smk3 > + +mdinc10c +mdinc10sq+ pwhtc +pwhtsq +edu2+ edu3 +husbgs > +husbcol+ husbmiss > + +currpmh +pastpmh +neverpmh, random= ~(1|id) > +(1|cogtest), data=global) > > Output indicating that this spline should probably look linear: > >> summary(global.gamm4$mer) > Random effects: > Groups Name Variance Std.Dev. > id (Intercept) 0.1823454 0.427019 > cogtest (Intercept) 0.0025498 0.050496 > Xr.1 s(adjtibx) 0.0000000 0.000000 > Residual 0.7782969 0.882211 > > Xs(adjtibx)Fx1 -0.0387360 0.0215596 -1.797 > > > Output getting a df for this spline of 0.20. > >> summary(global.gamm4$gam) > > Approximate significance of smooth terms: > edf Ref.df F p-value > s(adjtibx) 0.2009 0.2009 16.07 NA > > The plot looks linear, but reports a df =0.20. > > > So...to summarize my questions: > > 1. Are the splines produced by s(exp, fx=FALSE) or s(exp, fx=TRUE, k=k) > correct even though the reported degrees of freedom appears to be wrong? > > 2. Can I believe my plot? > > 3. How can I get the true df used when I use s(exp, fx=FALSE)? > > > > > > Thanks for any and all help you can provide! > > Melinda > > [[alternative HTML version deleted]] > > ______________________________________________ > 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. >-- Simon Wood, Mathematical Science, University of Bath BA2 7AY UK +44 (0)1225 386603 http://people.bath.ac.uk/sw283
Simon Wood
2011-Aug-04 09:30 UTC
[R] Incorrect degrees of freedom for splines using GAMM4?
Thanks for reporting this. It was a bug, now fixed in gamm4 0.1-3. Simon On 19/07/11 22:16, Melinda Power wrote:> Hello, > > I'm running mixed models in GAMM4 with 2 (non-nested) random intercepts and > I want to include a spline term for one of my exposure variables. However, > when I include a spline term, I always get reported degrees of freedom of > less than 1, even when I know that my spline is using more than 1 degree of > freedom. For example, here is the code for my model: > >> global.gamm4<-gamm4(zcog~s(adjpatx, fx=TRUE, k=5)+int234+cogagec+cogagesq > + > + + oldfran +newus +alc2 +alc3 +alc4 +alcmiss +smk2 +smk3 > + +mdinc10c +mdinc10sq+ pwhtc +pwhtsq +edu2+ edu3 +husbgs > +husbcol+ husbmiss > + +currpmh +pastpmh +neverpmh, random= ~(1|id) > +(1|cogtest), data=global) > > Using> summary(global.gamm4$mer), I get the following output for my spline > term, indicating that I use the expected 4 degrees of freedom. > > Xs(adjpatx)Fx1 0.1018943 0.1073225 0.949 > Xs(adjpatx)Fx2 -0.0708114 0.1123845 -0.630 > Xs(adjpatx)Fx3 0.7459511 0.6836413 1.091 > Xs(adjpatx)Fx4 -0.2062321 0.0923569 -2.233 > > However, when I use> summary(global.gamm4$gam). I get an estimate of > degrees of freedom that is not 4: > > Approximate significance of smooth terms: > edf Ref.df F p-value > s(adjpatx) 0.7588 0.7588 1.346 0.234 > > This degree of freedom = 0.76 also shows up on my plot. > > Ultimately, I would like to use a cubic regression penalized spline, > allowing R to choose the degrees of freedom for me using GCV. However, when > I use the correct code for this or variants of it using mgcv, I also get > degrees of freedom less than 1. For example, in the following code provides > a degree of freedom of less than 1 as well: > > >> global.gamm4<-gamm4(zcog~s(adjpatx, fx=FALSE)+int234+cogagec+cogagesq + > + + oldfran +newus +alc2 +alc3 +alc4 +alcmiss +smk2 +smk3 > + +mdinc10c +mdinc10sq+ pwhtc +pwhtsq +edu2+ edu3 +husbgs > +husbcol+ husbmiss > + +currpmh +pastpmh +neverpmh, random= ~(1|id) > +(1|cogtest), data=global) > > Output indicating that this spline should probably look linear: > >> summary(global.gamm4$mer) > Random effects: > Groups Name Variance Std.Dev. > id (Intercept) 0.1823454 0.427019 > cogtest (Intercept) 0.0025498 0.050496 > Xr.1 s(adjtibx) 0.0000000 0.000000 > Residual 0.7782969 0.882211 > > Xs(adjtibx)Fx1 -0.0387360 0.0215596 -1.797 > > > Output getting a df for this spline of 0.20. > >> summary(global.gamm4$gam) > > Approximate significance of smooth terms: > edf Ref.df F p-value > s(adjtibx) 0.2009 0.2009 16.07 NA > > The plot looks linear, but reports a df =0.20. > > > So...to summarize my questions: > > 1. Are the splines produced by s(exp, fx=FALSE) or s(exp, fx=TRUE, k=k) > correct even though the reported degrees of freedom appears to be wrong? > > 2. Can I believe my plot? > > 3. How can I get the true df used when I use s(exp, fx=FALSE)? > > > > > > Thanks for any and all help you can provide! > > Melinda > > [[alternative HTML version deleted]] > > ______________________________________________ > 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. >-- Simon Wood, Mathematical Science, University of Bath BA2 7AY UK +44 (0)1225 386603 http://people.bath.ac.uk/sw283