[R] What is the degrees of freedom in an nlme model

Tue Jul 13 17:20:14 CEST 2010

If the curves are sufficiently close to sine (cosine) curves and you know the period, then this can be restructured as a linear model and you can avoid all the complexities that come with non-linear models.  Further, from your description, it does not sound like you really gain much from using the mixed effects vs. just fixed effects, so this could be reduced to a simple use of lm and anova.

Hope this helps,

-- 
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
greg.snow at imail.org
801.408.8111

> -----Original Message-----
> From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-
> project.org] On Behalf Of Jun Shen
> Sent: Monday, July 12, 2010 4:57 PM
> To: Bert Gunter
> Cc: R-help
> Subject: Re: [R] What is the degrees of freedom in an nlme model
> 
> Hi, Bert,
> 
> Thanks for your thoughtful explanation. I think the problem is quite
> over my head and maybe I should leave "how" for experts :)
> 
> The situation is I have a group of sigmoid curves (let's say, they are
> supposed to be the same) but occasionally you will see a few curves
> kind of different. So how do we say they are actually different or not
> from the majority curves in a statistical way? The original idea was
> proposed by Monson and Rodbard in 1978 (Am J Physiol. 1978
> Aug;235(2):E97-102). The paper is freely available. The idea is to fit
> the curves individually and obtain the residual sum of squares and
> then fit the curves altogether somehow constraining some parameters
> and then you have another residual sum of squares. Then you can do a
> F-test. In my case, I wonder if I can use a mixed-effect modeling to
> do the simultaneous fitting job. Now you see, the problem is the
> degrees of freedom. Based on your explanation, it seems no reliable
> calculation of df for nonlinear models. However I can still see the df
> reported in nlme or nls models. Now I am not even sure if I should use
> them.
> 
> Another thing I observed is even I added more random effects to the
> nlme model, the denominator df did not seem to change. Is it correct?
> Thanks again.
> 
> Jun
> 
> On Mon, Jul 12, 2010 at 4:00 PM, Bert Gunter <gunter.berton at gene.com>
> wrote:
> > Jun:
> >
> > Short answer: There is no such thing as df for a nonlinear model
> (whether or
> > not mixed effects).
> >
> > Longer answer: df is the dimension of the null space when the data
> are
> > projected on the linear subspace of the model matrix of a **linear
> model **
> > . So, strictly speaking, no linear model, no df.
> >
> > HOWEVER... nonlinear models are usually (always??) fit by successive
> linear
> > approximations, and approximate df are obtained from these
> approximating
> > subspaces.
> >
> > However, the problem with this is that there is no guarantee that the
> > relevant residual distributions are sufficiently chisq with the
> approximate
> > df to give reasonable answers. In fact, lots of people much smarter
> than I
> > have spent lots of time trying to figure out what sorts of
> approximations
> > one should use to get trustworthy results. The thing is, in nonlinear
> > models, it can DEPEND on the exact form of the model -- indeed,
> that's what
> > distinguishes nonlinear models from linear ones! So this turns out to
> be
> > really hard and afaik these smart people don't agree on what should
> be done.
> >
> >
> > To see what one of the smartest people have to say about this, search
> the
> > archives for Doug Bates's comments on this w.r.t. lmer (he won't
> compute
> > such distributions nor provide P values because he doesn't know how
> to do it
> > reliably. Doug -- please correct me if I have it wrong).
> >
> > A stock way to extricate oneself from this dilemma is: bootstrap!
> > Unfortunately, this is also probably too facile: for one thing, with
> a
> > nondiagonal covariance matrix (as in mixed effects models), how do
> you
> > resample to preserve the covariance structure? I believe this is an
> area of
> > active research in the time series literature, for example. For
> another,
> > this may be too computationally demanding to be practicable due to
> > convergence issues.
> >
> > Bottom line: there may be no good way to do what you want.
> >
> > Note to experts: Please view this post as an invitation to correct my
> errors
> > and provide authoritative info.
> >
> > Cheers to all,
> >
> > Bert
> >
> > Bert Gunter
> > Genentech Nonclinical Biostatistics
> >
> >
> >
> > -----Original Message-----
> > From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-
> project.org] On
> > Behalf Of Jun Shen
> > Sent: Monday, July 12, 2010 12:34 PM
> > To: R-help
> > Subject: [R] What is the degrees of freedom in an nlme model
> >
> > Dear all,
> >
> > I want to do a F test, which involves calculation of the degrees of
> > freedom for the residuals. Now say, I have a nlme object "mod.nlme".
> I
> > have two questions
> >
> > 1.How do I extract the degrees of freedom?
> > 2.How is this degrees of freedom calculated in an nlme model?
> >
> > Thanks.
> >
> > Jun Shen
> >
> > Some sample code and data
> > =================================================================
> > mod.nlme<-nlme(RESP~E0+(Emax-
> E0)*CP**gamma/(EC50**gamma+CP**gamma),data=Data
> > ,
> > fixed=E0+Emax+gamma+EC50~1,
> > random=list(pdDiag(EC50+E0+gamma~1)),
> > groups=~ID,
> > start=list(fixed=c(E0=1,Emax=100,gamma=1,b=50))
> > )
> >
> > The "Data" object
> > structure(list(ID = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
> > 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
> > 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
> > 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
> > 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L,
> > 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L,
> > 6L, 6L, 6L, 6L, 6L, 6L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L,
> > 7L, 7L, 7L, 7L, 7L, 7L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L,
> > 8L, 8L, 8L), CP = c(1, 2, 3, 4.5, 5, 7.5, 11.25, 12, 18, 30,
> > 45, 60, 90, 120, 150, 200, 1, 2, 3, 4.5, 5, 7.5, 11.25, 12, 18,
> > 30, 45, 60, 90, 120, 150, 200, 1, 2, 3, 4.5, 5, 7.5, 11.25, 12,
> > 18, 30, 45, 60, 90, 120, 150, 200, 1, 2, 3, 4.5, 5, 7.5, 11.25,
> > 12, 18, 30, 45, 60, 90, 120, 150, 200, 1, 2, 3, 4.5, 5, 7.5,
> > 11.25, 12, 18, 30, 45, 60, 90, 120, 150, 200, 1, 2, 3, 4.5, 5,
> > 7.5, 11.25, 12, 18, 30, 45, 60, 90, 120, 150, 200, 1, 2, 3, 4.5,
> > 5, 7.5, 11.25, 12, 18, 30, 45, 60, 90, 120, 150, 200, 1, 2, 3,
> > 4.5, 5, 7.5, 11.25, 12, 18, 30, 45, 60, 90), RESP = c(3.19, 2.52,
> > 2.89, 3.28, 3.82, 7.15, 11.2, 16.25, 30.32, 55.25, 73.56, 82.07,
> > 89.08, 95.86, 97.97, 99.03, 3.49, 4.4, 3.54, 4.99, 3.81, 10.12,
> > 21.59, 24.93, 40.18, 61.01, 78.65, 88.81, 93.1, 94.61, 98.83,
> > 97.86, 0.42, 0, 2.58, 5.67, 3.64, 8.01, 12.75, 13.27, 24.65,
> > 46.1, 65.16, 77.74, 87.99, 94.4, 96.05, 100.4, 2.43, 0, 6.32,
> > 5.59, 8.48, 12.32, 26.4, 28.36, 43.38, 69.56, 82.53, 91.36, 95.37,
> > 98.36, 98.66, 98.8, 5.16, 2, 5.65, 3.48, 5.78, 5.5, 11.55, 8.53,
> > 18.02, 38.11, 58.93, 70.93, 85.62, 89.53, 96.19, 96.19, 2.76,
> > 2.99, 3.75, 3.02, 5.44, 3.08, 8.31, 10.85, 13.79, 32.06, 50.22,
> > 63.7, 81.34, 89.59, 93.06, 92.47, 3.32, 1.14, 2.43, 2.75, 3.02,
> > 5.4, 8.49, 7.91, 15.17, 35.01, 53.91, 68.51, 83.12, 86.85, 92.17,
> > 95.72, 3.58, 0.02, 3.69, 4.34, 6.32, 5.15, 9.7, 11.39, 23.38,
> > 42.9, 61.91, 71.82, 87.83)), .Names = c("ID", "CP", "RESP"), class =
> > "data.frame", row.names = c(NA,
> > -125L))
> >
> > ______________________________________________
> > 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.
> >
> >
> 
> ______________________________________________
> 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.