[R-sig-ME] two dependent variables
Iasonas Lamprianou
lamprianou at yahoo.com
Thu Nov 22 07:14:39 CET 2007
Dear friends, I have two binomial dependent variables. What sort of analysis can I use? Which package?
Jason
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: Thursday, 15 November, 2007 1:00:01 PM
Subject: R-sig-mixed-models Digest, Vol 11, Issue 8
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. Resid() and estimable() functions with lmer (Gustaf Granath)
2. Re: Resid() and estimable() functions with lmer (Douglas Bates)
3. non-linear GLMMs (Irene Mantzouni)
4. mixed model with crossed effects (Sharif S. Aly)
----------------------------------------------------------------------
Message: 1
Date: Wed, 14 Nov 2007 12:58:28 +0100
From: Gustaf Granath <Gustaf.Granath at ebc.uu.se>
Subject: [R-sig-ME] Resid() and estimable() functions with lmer
To: r-sig-mixed-models at r-project.org
Message-ID: <473AE2E4.9030205 at ebc.uu.se>
Content-Type: text/plain; charset=ISO-8859-1; format=flowed
Hi all,
Two questions:
1. Is there a way to evaluate models from lmer() with a poisson
distribution? I get the following error message:
library(lme4)
lmer(tot.fruit~infl.treat+def.treat+(1|initial.size),family=poisson)->model
plot(fitted(model),resid(model))
Error: 'resid' is not implemented yet
Are there any other options?
2. Why doesn't the function estimable() in gmodels work with lmer()
using a poisson distribution? I know that my coefficients are right
because it works fine if I run the model without poisson distribution.
#CODE (with latest R version and package updates)#
#Without Poisson distribution#
library(lme4)
library(gmodels)
lmer(tot.fruit~infl.treat+def.treat+(1|initial.size))->model1
summary(model1)
Linear mixed-effects model fit by REML
Formula: tot.fruit ~ infl.treat + def.treat + (1 | initial.size)
AIC BIC logLik MLdeviance REMLdeviance
2449 2471 -1218 2447 2437
Random effects:
Groups Name Variance Std.Dev.
initial.size (Intercept) 71.585 8.4608
Residual 49.856 7.0608
number of obs: 323, groups: initial.size, 292
Fixed effects:
Estimate Std. Error t value
(Intercept) 12.8846 1.3028 9.890
infl.treat1 -0.4738 1.1819 -0.401
def.treat2 -3.5522 1.6022 -2.217
def.treat3 -2.1757 1.6461 -1.322
def.treat4 -2.1613 1.7003 -1.271
Correlation of Fixed Effects:
(Intr) infl.1 df.tr2 df.tr3
infl.treat1 -0.413
def.treat2 -0.616 0.013
def.treat3 -0.641 0.002 0.493
def.treat4 -0.638 0.028 0.469 0.524
#Coefficients#
mean.infl.treat0=c(1,0,1/4,1/4,1/4)
mean.infl.treat1=c(1,1,1/4,1/4,1/4)
mean.def.treat1=c(1,1/2,0,0,0)
mean.def.treat2=c(1,1/2,1,0,0)
mean.def.treat3=c(1,1/2,0,1,0)
mean.def.treat4=c(1,1/2,0,0,1)
means=rbind(mean.infl.treat0,mean.infl.treat1,mean.def.treat1,mean.def.treat2,mean.def.treat3,mean.def.treat4)
estimable(model1,means,sim.lmer=TRUE,n.sim=1000)
Estimate Std. Error p value
(1 0 0.25 0.25 0.25) 10.897825 0.8356260 0
(1 1 0.25 0.25 0.25) 10.421633 0.9364839 0
(1 0.5 0 0 0) 12.577408 1.2055979 0
(1 0.5 1 0 0) 9.097063 1.1769760 0
(1 0.5 0 1 0) 10.523485 1.2238635 0
(1 0.5 0 0 1) 10.440959 1.2031459 0
#With Poisson distribution#
lmer(tot.fruit~infl.treat+def.treat+(1|initial.size),family=poisson)->model2
summary(model2)
Generalized linear mixed model fit using Laplace
Formula: tot.fruit ~ infl.treat + def.treat + (1 | initial.size)
Family: poisson(log link)
AIC BIC logLik deviance
1073 1096 -530.5 1061
Random effects:
Groups Name Variance Std.Dev.
initial.size (Intercept) 0.84201 0.91761
number of obs: 323, groups: initial.size, 292
Estimated scale (compare to 1 ) 1.130529
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.08865 0.10478 19.934 < 2e-16 ***
infl.treat1 0.10615 0.09412 1.128 0.25941
def.treat2 -0.37354 0.11398 -3.277 0.00105 **
def.treat3 -0.18566 0.12679 -1.464 0.14310
def.treat4 -0.10624 0.13451 -0.790 0.42962
---
Signif. codes: 0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1
Correlation of Fixed Effects:
(Intr) infl.1 df.tr2 df.tr3
infl.treat1 -0.418
def.treat2 -0.554 0.016
def.treat3 -0.599 -0.001 0.467
def.treat4 -0.623 0.055 0.460 0.548
estimable(model2,lsmeans,sim.lmer=TRUE,n.sim=1000)
Error in FUN(newX[, i], ...) :
`param' has no names and does not match number of coefficients of
model. Unable to construct coefficient vector
#End of CODE#
Any ideas as to why this is? It does say in the help of ?estimable
that it should work for lmer objects.
Thanks in advance for any help,
Gustaf (PhD student)
Dept of Plant Ecology
Evolutionary Biology Centre (EBC)
Uppsala University
------------------------------
Message: 2
Date: Wed, 14 Nov 2007 07:17:45 -0600
From: "Douglas Bates" <bates at stat.wisc.edu>
Subject: Re: [R-sig-ME] Resid() and estimable() functions with lmer
To: Gustaf.Granath at ebc.uu.se
Cc: r-sig-mixed-models at r-project.org
Message-ID:
<40e66e0b0711140517l69a470cx53929c094002144d at mail.gmail.com>
Content-Type: text/plain; charset=ISO-8859-1
On Nov 14, 2007 5:58 AM, Gustaf Granath <Gustaf.Granath at ebc.uu.se> wrote:
> Hi all,
> Two questions:
> 1. Is there a way to evaluate models from lmer() with a poisson
> distribution? I get the following error message:
> library(lme4)
> lmer(tot.fruit~infl.treat+def.treat+(1|initial.size),family=poisson)->model
> plot(fitted(model),resid(model))
> Error: 'resid' is not implemented yet
The model can be fit, it is the resid function that is not yet
implemented. The reason it is not yet implemented is because there
are many different kinds of residuals for generalized linear mixed
models. It requires considerably more design and coding than does
obtaining the residuals for a linear mixed model, and even those are
difficult to define for the general model (Do you incorporate the
random effects or not? If you do, what do you use in the case of
multiple grouping factors?)
Right now I am thinking about other issues in the form of the model so
I can incorporate both generalized linear mixed models and nonlinear
mixed models. Residuals for generalized linear mixed models are
further down the priority list.
> Are there any other options?
>
> 2. Why doesn't the function estimable() in gmodels work with lmer()
> using a poisson distribution? I know that my coefficients are right
> because it works fine if I run the model without poisson distribution.
I don't know what the estimable function expects the structure of an
lmer or glmer model to be but it would not surprise me if it did not
coincide with the current version. The structure is continuously
evolving, not because I want to make things difficult for others but
because my understanding of the model and the computational methods
evolves. On something as complicated as a generalized nonlinear mixed
model with crossed grouping factors it is not surprising that one
doesn't get the design down on the first try. At least this one
doesn't. I hope the structure of the classes and the algorithms will
stabilize soon (say, by the end of the year).
Even then it will be difficult to implement something like estimable
without doing finite difference calculations for generalized linear
mixed models. A matrix that determines the marginal variance
covariance for the fixed-effects estimates (well, not quite but I
don't want to get into the technicalities) is calculated for linear
mixed models but not for generalized linear mixed models.
> #CODE (with latest R version and package updates)#
> #Without Poisson distribution#
>
> library(lme4)
> library(gmodels)
> lmer(tot.fruit~infl.treat+def.treat+(1|initial.size))->model1
> summary(model1)
>
> Linear mixed-effects model fit by REML
> Formula: tot.fruit ~ infl.treat + def.treat + (1 | initial.size)
> AIC BIC logLik MLdeviance REMLdeviance
> 2449 2471 -1218 2447 2437
> Random effects:
> Groups Name Variance Std.Dev.
> initial.size (Intercept) 71.585 8.4608
> Residual 49.856 7.0608
> number of obs: 323, groups: initial.size, 292
>
> Fixed effects:
> Estimate Std. Error t value
> (Intercept) 12.8846 1.3028 9.890
> infl.treat1 -0.4738 1.1819 -0.401
> def.treat2 -3.5522 1.6022 -2.217
> def.treat3 -2.1757 1.6461 -1.322
> def.treat4 -2.1613 1.7003 -1.271
>
> Correlation of Fixed Effects:
> (Intr) infl.1 df.tr2 df.tr3
> infl.treat1 -0.413
> def.treat2 -0.616 0.013
> def.treat3 -0.641 0.002 0.493
> def.treat4 -0.638 0.028 0.469 0.524
>
> #Coefficients#
>
> mean.infl.treat0=c(1,0,1/4,1/4,1/4)
> mean.infl.treat1=c(1,1,1/4,1/4,1/4)
> mean.def.treat1=c(1,1/2,0,0,0)
> mean.def.treat2=c(1,1/2,1,0,0)
> mean.def.treat3=c(1,1/2,0,1,0)
> mean.def.treat4=c(1,1/2,0,0,1)
> means=rbind(mean.infl.treat0,mean.infl.treat1,mean.def.treat1,mean.def.treat2,mean.def.treat3,mean.def.treat4)
> estimable(model1,means,sim.lmer=TRUE,n.sim=1000)
>
> Estimate Std. Error p value
> (1 0 0.25 0.25 0.25) 10.897825 0.8356260 0
> (1 1 0.25 0.25 0.25) 10.421633 0.9364839 0
> (1 0.5 0 0 0) 12.577408 1.2055979 0
> (1 0.5 1 0 0) 9.097063 1.1769760 0
> (1 0.5 0 1 0) 10.523485 1.2238635 0
> (1 0.5 0 0 1) 10.440959 1.2031459 0
>
> #With Poisson distribution#
>
> lmer(tot.fruit~infl.treat+def.treat+(1|initial.size),family=poisson)->model2
> summary(model2)
>
> Generalized linear mixed model fit using Laplace
> Formula: tot.fruit ~ infl.treat + def.treat + (1 | initial.size)
> Family: poisson(log link)
> AIC BIC logLik deviance
> 1073 1096 -530.5 1061
> Random effects:
> Groups Name Variance Std.Dev.
> initial.size (Intercept) 0.84201 0.91761
> number of obs: 323, groups: initial.size, 292
>
> Estimated scale (compare to 1 ) 1.130529
>
> Fixed effects:
> Estimate Std. Error z value Pr(>|z|)
> (Intercept) 2.08865 0.10478 19.934 < 2e-16 ***
> infl.treat1 0.10615 0.09412 1.128 0.25941
> def.treat2 -0.37354 0.11398 -3.277 0.00105 **
> def.treat3 -0.18566 0.12679 -1.464 0.14310
> def.treat4 -0.10624 0.13451 -0.790 0.42962
> ---
> Signif. codes: 0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1
>
> Correlation of Fixed Effects:
> (Intr) infl.1 df.tr2 df.tr3
> infl.treat1 -0.418
> def.treat2 -0.554 0.016
> def.treat3 -0.599 -0.001 0.467
> def.treat4 -0.623 0.055 0.460 0.548
>
> estimable(model2,lsmeans,sim.lmer=TRUE,n.sim=1000)
> Error in FUN(newX[, i], ...) :
> `param' has no names and does not match number of coefficients of
> model. Unable to construct coefficient vector
>
> #End of CODE#
>
> Any ideas as to why this is? It does say in the help of ?estimable
> that it should work for lmer objects.
>
> Thanks in advance for any help,
>
> Gustaf (PhD student)
>
>
> Dept of Plant Ecology
> Evolutionary Biology Centre (EBC)
> Uppsala University
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
------------------------------
Message: 3
Date: Wed, 14 Nov 2007 20:23:52 +0100
From: "Irene Mantzouni" <ima at difres.dk>
Subject: [R-sig-ME] non-linear GLMMs
To: <R-SIG-Mixed-Models at r-project.org>
Message-ID:
<68E7981938EAF54F987AD3848A0A6416E5838B at ka-mail01.dfu.local>
Content-Type: text/plain; charset="ISO-8859-7"
Is it (or will it become...) possible to fit a mixed non-linear SSmicmen model assuming log-normal error?
I have tried to use a log transformed micmen model but it does not seem so powerful.
Any other ideas about possible approaches are more than welcome!
Irene
------------------------------
Message: 4
Date: Wed, 14 Nov 2007 23:02:40 -0800
From: "Sharif S. Aly" <saly at ucdavis.edu>
Subject: [R-sig-ME] mixed model with crossed effects
To: r-sig-mixed-models at r-project.org
Message-ID: <473BEF10.3070105 at ucdavis.edu>
Content-Type: text/plain; charset=ISO-8859-1; format=flowed
Dear list,
If you please help me with the R command for a 4 level random effects
model I would be grateful. I chose to put the details of my study at
the end of the email as some may prefer not to read it.
My outcome is continuous (normally distributed). All levels are factor
levels. The lowest level (veterinarian) is crossed with level 2 (day of
sampling) which is crossed with level 3 (sampled pen or location on
farm) which is nested in level 4 (dairy).
I understand I have to specify the grouping structure of the data before
the model command, is that correct? My model so far is:
XCoutcome<-GroupedData: outcome~ Dairy | Pen.number | cons
lme(AvgPCR~ 1, random = ~ 1 | Dairy /Pen.number /
pdBlocked(list(pdIdent(~ENV),pdIdent(~Veterinarian))), data=XCoutcome)
Best regards,
Sharif
Description of study design:
My data is comprised of bacteria counts in samples collected by 2
veterinarians over 3 days from different pens (stalls that house cows)
nested in different dairies. Vet is crossed with day which is crossed
with pen (because vet 1 is the same vet 1 who sampled on all days (as is
vet 2), and day 1 is the same day 1 of sampling in all pens (as is days
2 and 3)); pens however are nested in dairies, meaning that pen 1 in
dairy 1 is different that pen 1 in dairy 2, 3 and 4.
My objective was to estimate the similarity in bacterial counts in
samples collected by 2 vets, in the same pen, in the same dairy (that
specific intraclass correlation coefficient) and for reasons I can
explain, I chose to have only a fixed effect intercept (basically we are
not interested in the effect of any particular vet, pen or dairy or day,
i.e. we do not wish to estimate the effect of Nov 16 per say)
Sharif Aly,
Graduate group in Epidemiology,
University of California, Davis
------------------------------
_______________________________________________
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 11, Issue 8
*************************************************
___________________________________________________________
Yahoo! Answers - Got a question? Someone out there knows the answer. Try it
now.
http://uk.answers.yahoo.com/
More information about the R-sig-mixed-models
mailing list