[R-sig-ME] R-sig-mixed-models
Iasonas Lamprianou
lamprianou at yahoo.com
Fri Oct 5 12:22:19 CEST 2007
Dear friends,
I am trying the following commands, because I want to explore my data before doing an analysis with TEACHER as random effect.
> data <-groupedData( DependentVariable ~ Predictor |SCHOOL/TEACHER,data=data)
> samp <- sample(levels(data$TEACHER),29)
> level2.subgroup <- subset(data, TEACHER %in% samp)
> level2<-lmList(DependentVariable ~ Predictor | TEACHER,data=level2.subgroup)
and i get this error Error in na.fail.default(data) : missing values in object
what is this error? What do I do wrong?
Dr. Iasonas Lamprianou
Department of Education
The University of Manchester
Oxford Road, Manchester M13 9PL, UK
Tel. 0044 161 275 3485
iasonas.lamprianou at manchester.ac.uk
----- Original Message ----
From: "r-sig-mixed-models-request at r-project.org" <r-sig-mixed-models-request at r-project.org>
To: r-sig-mixed-models at r-project.org
Sent: Friday, 5 October, 2007 1:41:25 AM
Subject: R-sig-mixed-models Digest, Vol 10, Issue 4
Send R-sig-mixed-models mailing list submissions to
r-sig-mixed-models at r-project.org
To subscribe or unsubscribe via the World Wide Web, visit
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
or, via email, send a message with subject or body 'help' to
r-sig-mixed-models-request at r-project.org
You can reach the person managing the list at
r-sig-mixed-models-owner at r-project.org
When replying, please edit your Subject line so it is more specific
than "Re: Contents of R-sig-mixed-models digest..."
Today's Topics:
1. Re: Gastric emptying data: nlmer vs. nlme (Douglas Bates)
2. lmer vs lmer2 (dave fournier)
3. Re: lmer vs lmer2 (Douglas Bates)
----------------------------------------------------------------------
Message: 1
Date: Thu, 4 Oct 2007 14:52:28 -0500
From: "Douglas Bates" <bates at stat.wisc.edu>
Subject: Re: [R-sig-ME] Gastric emptying data: nlmer vs. nlme
To: "Dieter Menne" <dieter.menne at menne-biomed.de>
Cc: r-sig-mixed-models at r-project.org
Message-ID:
<40e66e0b0710041252w73c73e4aj340e031d2db416a4 at mail.gmail.com>
Content-Type: text/plain; charset=ISO-8859-1
Dieter,
Thanks for the very thorough description of the problem and for
including a reproducible example. I enclose a modified version of
your script and the output using the development version of the lme4
package (https://svn.R-project.org/R-packages/branches/gappy-lmer/)
You will see that the development version does better on the first
model but still ends up giving a "false convergence" message. I then
fit ge0A.nlmer and added another model ge0B.nlmer that removes the
random effect for kappa. The (conservative) p-value for a test of H0:
ge0B.nlmer versus Ha: ge0A.nlmer is
> pchisq(979.76761 - 973.73274, df = 2)
[1] 0.9510734
so incorporating random effects for kappa is, at best, marginally significant.
The reason that the first model is so difficult to fit is because
there are 6 variance-covariance parameters and only 8 levels of
subj:treat. You are trying to estimate too many variance-covariance
parameters from too few groups. The likelihood surface will be very
flat and the parameter estimates will be ill-defined.
The reason that nlmer from the 0.99875-8 release of lme4 gave up has
to do with the calculation of the conditional modes of the random
effects to evaluate the Laplace approximation to the deviance. In
that version I retained the values of the conditional modes of the
random effects (these are the values the maximum the conditional
density of the random effects given the data and the current values of
the model parameters - the Laplace approximation is evaluated at these
values and the adaptive Gauss-Hermite approximation is centered around
these values) between evaluations of the deviance. That's a good idea
because the vector of the conditional modes for a different set of
parameters is going to be similar to the current values so you have
good starting estimates. However, it is a bad idea in that the
evaluation of the approximation to the deviance will not only depend
on the values of the parameters and the data but also on where the
last value was taken. This means that the function value being
optimized is not reproducible and that causes a lot of problems in a
derivative-free optimization.
To avoid this I now start each evaluation of the conditional modes at
the same point (all random effects start at zero) so the evaluation is
reproducible.
On 10/1/07, Dieter Menne <dieter.menne at menne-biomed.de> wrote:
> Dear Group,
>
> I have tried to redo published nlme-fits of gastric emptying data recorded by
> MRI with nlmer. nlme needs time, but give reasonable results, informing me about
> the known correlation between two nonlinear parameters, while nlmer gave up
> after two iterations without telling me why.
>
> I hope I got the translation to nlmer syntax correctly. What's wrong?
>
> Dieter Menne
>
>
> ----------------------------------------------------
> # Gastric Emptying data from
> # Goetze et. al (University Hospital of Z?rich)
> # http://ajpgi.physiology.org/cgi/content/abstract/00498.2005v1
> # Demo data set from http://www.menne-biomed.de/gastempt
> # R-Version: 2.5.1, i386, mingw32
> # lme4_0.99875-8.zip
>
> useNlme = TRUE
>
> # Get data
> ge = read.table("http://www.menne-biomed.de/gastempt/gastempt.csv";,
> sep=",",header=TRUE)
> #ge = read.table("gastempt.csv",sep=",",header=TRUE)
> ge$subj = as.factor(ge$subj)
> ge$subjtreat = as.factor(paste(ge$subj,ge$treat,sep="."))
> ge$t = as.double(ge$t)
>
>
> EmptInit= function(mCall,LHS,data){ # dummy, not explicitely used
> stop("Should not be called")
> }
>
> # Standard LinExp model for gastric emptying
>
> SSEmptLinExp=selfStart(~v0*(1+kappa*t/tempt)*exp(-t/tempt),
> initial=EmptInit, parameters= c("v0","kappa","tempt"))
>
> # nlme final is 643,1.4,74
> start = list(fixed=c(v0=600,kappa=1.0,tempt=70))
>
> if (useNlme) {
> library(nlme)
> contr = nlmeControl(pnlsTol=0.3,pnlsMaxIter=200,msMaxIter=200)
> ge0.nlme=nlme(v~SSEmptLinExp(t,v0,kappa,tempt),
> fixed =list(v0+kappa+tempt~1),
> groups = ~subjtreat,
> control=contr, start=start, data=ge, verbose=F
> )
> summary(ge0.nlme)
> # Note correlation between kappa and tempt
> ge0A.nlme = update(ge0.nlme,random=v0+kappa~1)
> summary(ge0A.nlme)
> anova(ge0A.nlme,ge0.nlme)
> }
>
> if (!useNlme) {
> library(lme4) # lme4_0.99875-8.zip
> # This stops after two iterations without telling me details
> ge0.nlmer = nlmer(
> v~SSEmptLinExp(t,v0,kappa,tempt)~(v0+kappa+tempt|subjtreat),
> data=ge, start=start,verbose=T)
> show(ge0.nlmer)
> # This is close to the result of ge0A.nlme)
> ge0A.nlmer = nlmer(
> v~SSEmptLinExp(t,v0,kappa,tempt)~(v0+kappa|subjtreat),
> data=ge, start=start,verbose=T)
> show(ge0.nlmer)
> }
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
------------------------------
Message: 2
Date: Thu, 04 Oct 2007 22:46:50 -0700
From: dave fournier <otter at otter-rsch.com>
Subject: [R-sig-ME] lmer vs lmer2
To: r-sig-mixed-models at r-project.org
Message-ID: <4705CFCA.9030100 at otter-rsch.com>
Content-Type: text/plain; charset=ISO-8859-1; format=flowed
Hi,
I checked this example out with ADMB-RE using a modification of
our glmmADMB program and have found the following:
1)
Parameter estimates with ADMB-RE are stable and
I get almost the same ones with or without the group 177 observations.
2) I get almost exactly the same LL estimate as SAS.
3) My estimates for the fixed effects are similar to those in
lmer2 except for the Intercept
Here are the estimates for lmer2 without group 177
Estimate Std. Error t value
(Intercept) -1.948119 0.095877 -20.32
Height 1.640650 0.032800 50.02
Age 0.019379 0.001310 14.79
InitHeight 0.143977 0.111043 1.30
InitAge -0.014618 0.007501 -1.95
these are the ADMB-RE estimates without group 177
LL = 2294.85
real_b -2.0369e+000 1.0393e-001
real_b 1.6460e+000 3.4587e-002
real_b 1.9275e-002 1.3685e-003
real_b 2.4857e-001 1.1984e-001
real_b -2.1290e-002 8.1749e-003
these are the estimates with group 177
real_b -2.0353e+000 1.0380e-001
real_b 1.6438e+000 3.4430e-002
real_b 1.9337e-002 1.3595e-003
real_b 2.5070e-001 1.1966e-001
real_b -2.1486e-002 8.1618e-003
Here are the lmer2 estimates with group 177 included
(Intercept) -2.048023 0.101413 -20.19
Height 1.643644 0.031106 52.84
Age 0.019092 0.001391 13.73
InitHeight 0.262909 0.118516 2.22
InitAge -0.021540 0.008111 -2.66
I think it is highly unlikely that the lmer2 estimate of
-1.948119 is the "correct" one and changes so much with
the addition of these few observations, while just by chance
ADMB-RE is wrong but happens to get the same estimate
for Intercept with and without group 177.
So it appears that lmer2 is not trustworthy.
Does anyone understand why the SAS point estimates appear to be
completely different?
Cheers,
Dave
David A. Fournier
P.O. Box 2040,
Sidney, B.C. V8l 3S3
Canada
Phone/FAX 250-655-3364
http://otter-
--
David A. Fournier
P.O. Box 2040,
Sidney, B.C. V8l 3S3
Canada
Phone/FAX 250-655-3364
http://otter-rsch.com
------------------------------
Message: 3
Date: Thu, 4 Oct 2007 17:41:11 -0500
From: "Douglas Bates" <bates at stat.wisc.edu>
Subject: Re: [R-sig-ME] lmer vs lmer2
To: davef at otter-rsch.com, "Garrett Fitzmaurice"
<fitzmaur at hsph.harvard.edu>, "Nan Laird" <LAIRD at hsph.harvard.edu>
Cc: r-sig-mixed-models at r-project.org
Message-ID:
<40e66e0b0710041541n486a8cbu687d1142b7062e94 at mail.gmail.com>
Content-Type: text/plain; charset="iso-8859-1"
On 10/5/07, dave fournier <otter at otter-rsch.com> wrote:
>
> Hi,
>
> I checked this example out with ADMB-RE using a modification of
> our glmmADMB program and have found the following:
>
> 1)
>
> Parameter estimates with ADMB-RE are stable and
> I get almost the same ones with or without the group 177 observations.
>
> 2) I get almost exactly the same LL estimate as SAS.
>
> 3) My estimates for the fixed effects are similar to those in
> lmer2 except for the Intercept
>
> Here are the estimates for lmer2 without group 177
> Estimate Std. Error t value
> (Intercept) -1.948119 0.095877 -20.32
> Height 1.640650 0.032800 50.02
> Age 0.019379 0.001310 14.79
> InitHeight 0.143977 0.111043 1.30
> InitAge -0.014618 0.007501 -1.95
>
> these are the ADMB-RE estimates without group 177
> LL = 2294.85
> real_b -2.0369e+000 1.0393e-001
> real_b 1.6460e+000 3.4587e-002
> real_b 1.9275e-002 1.3685e-003
> real_b 2.4857e-001 1.1984e-001
> real_b -2.1290e-002 8.1749e-003
>
> these are the estimates with group 177
>
> real_b -2.0353e+000 1.0380e-001
> real_b 1.6438e+000 3.4430e-002
> real_b 1.9337e-002 1.3595e-003
> real_b 2.5070e-001 1.1966e-001
> real_b -2.1486e-002 8.1618e-003
>
> Here are the lmer2 estimates with group 177 included
> (Intercept) -2.048023 0.101413 -20.19
> Height 1.643644 0.031106 52.84
> Age 0.019092 0.001391 13.73
> InitHeight 0.262909 0.118516 2.22
> InitAge -0.021540 0.008111 -2.66
>
> I think it is highly unlikely that the lmer2 estimate of
> -1.948119 is the "correct" one and changes so much with
> the addition of these few observations, while just by chance
> ADMB-RE is wrong but happens to get the same estimate
> for Intercept with and without group 177.
> So it appears that lmer2 is not trustworthy.
>
> Does anyone understand why the SAS point estimates appear to be
> completely different?
Because the SAS program is fitting a different model?
If you look at the sample SAS programs on the web site for the book
you will see that the authors are fitting models with fixed effects
for the logarithm of the height and the logarithm of the base height.
I have sort of lost track of the discussion of this example but I can
reproduce the results from Garrett Fitzmaurice's SAS analysis of these
data except for the variance-covariance of the random effects in the
model with correlated random effects for the intercept, the age and
the logarithm of the height. With the development version of the lme4
package I get a (near) singular variance-covariance matrix in that
model fit while SAS PROC MIXED doesn't indicate a problem with the
fit. The only indication of a problem from SAS is the large standard
errors on the estimates of the variance-covariance parameters.
I enclose the R script and output using the development version of the
lme4 package. I have copied the variable names, etc. from the SAS
programs on Garrett's web site. I fit two versions of each model, one
with all the subjects' data (fm1, fm2 and fm3) and one eliminating the
data for subject 197 (fm1a, fm2a and fm3a). (Dave: according to the
information on Garrett's web site it is subject 197, not 177, who
appears to be an outlier.)
The clue that model fm3a has a singular variance covariance matrix is
the estimated correlation of -1.000. Also, the verbose output shows
the converged value of the second parameter is very close to zero.
The first three parameters represent the variances of linear
combinations of the random effects. The interpretation is that a
linear combination of the random effects for the intercept and for age
has zero variance.
The big change in the development version of the lme4 package relative
to earlier versions is a rewriting of the mixed model equations so
that a singular variance-covariance matrix for the random effects is
approached smoothly, even though it is on the boundary.
I have permission from the book's authors to create an R package with
the data sets from the book. The package will be called AppLong and
will include sample analyses reproducing the SAS analyses as best I
can.
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: fev1_Rout.txt
Url: https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20071004/fd58543a/attachment.txt
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: fev1_R.txt
Url: https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20071004/fd58543a/attachment-0001.txt
------------------------------
_______________________________________________
R-sig-mixed-models mailing list
R-sig-mixed-models at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
End of R-sig-mixed-models Digest, Vol 10, Issue 4
*************************************************
___________________________________________________________
Want ideas for reducing your carbon footprint? Visit Yahoo! For Good http://uk.promotions.yahoo.com/forgood/environment.html
More information about the R-sig-mixed-models
mailing list