Hi, I run the following models: 1a. lmer(Y~X+(1|Subject),family=binomial(link="logit")) and 1b. lmer(Y~X+(1|Subject),family=binomial(link="logit"),method="PQL") Why does 1b produce results different from 1a? The reason why I am asking is that the help states that "PQL" is the default of GLMMs and 2. gamm(Y~X,family=binomial(link="logit"),random=list(Subject=~1)) The interesting thing about the example below is, that gamm is also supposed to fit by "PQL". Interestingly, however, the GAMM fit yields about the coefficient estimates of 1b. But the significance values of 1a. Any insight would be greatly appreciated. library(lme4) library(mgcv) Y=c(0,1,1,1,1,0,0,0,0,0,1,1,1,1,0,0,0,1,1,1,1) X=c(1,2,3,4,3,1,0,0,2,3,3,2,4,3,2,1,1,3,4,2,3) Subject=as.factor(c(1,2,3,4,5,6,7,1,2,3,4,5,6,7,1,2,3,4,5,6,7)) cbind(Y,X,Subject) r1=lmer(Y~X+(1|Subject),family=binomial(link="logit")) summary(r1) r2=lmer(Y~X+(1|Subject),family=binomial(link="logit"),method="PQL") summary(r2) r3=gamm(Y~X,family=binomial(link="logit"),random=list(Subject=~1)) summary(r3$gam) ------------------------- cuncta stricte discussurus
your r2 model corresponds to method = "Laplace". r4=lmer(Y~X+(1|Subject),family=binomial(link="logit"),method="Laplace") is equivalent to r2. Bests, Abderrahim ----- Original Message ----- From: "Daniel Malter" <daniel at umd.edu> To: <r-help at stat.math.ethz.ch> Sent: Friday, February 15, 2008 12:50 AM Subject: [R] LMER> Hi, > > I run the following models: > > 1a. lmer(Y~X+(1|Subject),family=binomial(link="logit")) and > 1b. lmer(Y~X+(1|Subject),family=binomial(link="logit"),method="PQL") > > Why does 1b produce results different from 1a? The reason why I am asking > is > that the help states that "PQL" is the default of GLMMs > > and > > 2. gamm(Y~X,family=binomial(link="logit"),random=list(Subject=~1)) > > The interesting thing about the example below is, that gamm is also > supposed > to fit by "PQL". Interestingly, however, the GAMM fit yields about the > coefficient estimates of 1b. But the significance values of 1a. Any > insight > would be greatly appreciated. > > > library(lme4) > library(mgcv) > > Y=c(0,1,1,1,1,0,0,0,0,0,1,1,1,1,0,0,0,1,1,1,1) > X=c(1,2,3,4,3,1,0,0,2,3,3,2,4,3,2,1,1,3,4,2,3) > Subject=as.factor(c(1,2,3,4,5,6,7,1,2,3,4,5,6,7,1,2,3,4,5,6,7)) > cbind(Y,X,Subject) > > r1=lmer(Y~X+(1|Subject),family=binomial(link="logit")) > summary(r1) > > r2=lmer(Y~X+(1|Subject),family=binomial(link="logit"),method="PQL") > summary(r2) > > r3=gamm(Y~X,family=binomial(link="logit"),random=list(Subject=~1)) > summary(r3$gam) > > > > ------------------------- > cuncta stricte discussurus > > ______________________________________________ > 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. >
Could you send us the output of sessionInfo() please so we can see
which version of the lme4 package you are using? In recent versions,
especially the development version available as
install.packages("lme4", repos =
"http://r-forge.r-project.org")
the PQL algorithm is no longer used. The Laplace approximation is now
the default. The adaptive Gauss-Hermite quadrature (AGQ)
approximation may be offered in the future.
If the documentation indicates that PQL is the default then that is a
documentation error. With the currently available implementation of
the direct optimization of the Laplace approximation to the
log-likelihood for the model there is no purpose in offering PQL.
On Thu, Feb 14, 2008 at 6:50 PM, Daniel Malter <daniel at umd.edu>
wrote:> Hi,
>
> I run the following models:
>
> 1a. lmer(Y~X+(1|Subject),family=binomial(link="logit")) and
> 1b.
lmer(Y~X+(1|Subject),family=binomial(link="logit"),method="PQL")
>
> Why does 1b produce results different from 1a? The reason why I am asking
is
> that the help states that "PQL" is the default of GLMMs
>
> and
>
> 2.
gamm(Y~X,family=binomial(link="logit"),random=list(Subject=~1))
>
> The interesting thing about the example below is, that gamm is also
supposed
> to fit by "PQL". Interestingly, however, the GAMM fit yields
about the
> coefficient estimates of 1b. But the significance values of 1a. Any
insight
> would be greatly appreciated.
>
>
> library(lme4)
> library(mgcv)
>
> Y=c(0,1,1,1,1,0,0,0,0,0,1,1,1,1,0,0,0,1,1,1,1)
> X=c(1,2,3,4,3,1,0,0,2,3,3,2,4,3,2,1,1,3,4,2,3)
> Subject=as.factor(c(1,2,3,4,5,6,7,1,2,3,4,5,6,7,1,2,3,4,5,6,7))
> cbind(Y,X,Subject)
>
> r1=lmer(Y~X+(1|Subject),family=binomial(link="logit"))
> summary(r1)
>
>
r2=lmer(Y~X+(1|Subject),family=binomial(link="logit"),method="PQL")
> summary(r2)
>
>
r3=gamm(Y~X,family=binomial(link="logit"),random=list(Subject=~1))
> summary(r3$gam)
>
>
>
> -------------------------
> cuncta stricte discussurus
>
> ______________________________________________
> 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.
>
Thanks for your replies. My real problem is that, for my real data, I get
basically the same results from r2 and r3 (so to speak), but the coefficient
estimates and significance levels for r1 are very different from those of r2
and r3. And therefore, I do not know which of the results to trust and which
not (if any).
The session info follows:
R version 2.6.0 (2007-10-03)
i386-pc-mingw32
locale:
LC_COLLATE=English_United States.1252;LC_CTYPE=English_United
States.1252;LC_MONETARY=English_United
States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] nlme_3.1-86 mgcv_1.3-29 lme4_0.99875-9 Matrix_0.999375-3
lattice_0.16-5
loaded via a namespace (and not attached):
[1] grid_2.6.0
Cheers,
Daniel
-------------------------
cuncta stricte discussurus
-------------------------
-----Urspr?ngliche Nachricht-----
Von: dmbates at gmail.com [mailto:dmbates at gmail.com] Im Auftrag von Douglas
Bates
Gesendet: Friday, February 15, 2008 7:29 AM
An: Daniel Malter
Cc: r-help at stat.math.ethz.ch
Betreff: Re: [R] LMER
Could you send us the output of sessionInfo() please so we can see which
version of the lme4 package you are using? In recent versions, especially
the development version available as
install.packages("lme4", repos =
"http://r-forge.r-project.org")
the PQL algorithm is no longer used. The Laplace approximation is now the
default. The adaptive Gauss-Hermite quadrature (AGQ) approximation may be
offered in the future.
If the documentation indicates that PQL is the default then that is a
documentation error. With the currently available implementation of the
direct optimization of the Laplace approximation to the log-likelihood for
the model there is no purpose in offering PQL.
On Thu, Feb 14, 2008 at 6:50 PM, Daniel Malter <daniel at umd.edu>
wrote:> Hi,
>
> I run the following models:
>
> 1a. lmer(Y~X+(1|Subject),family=binomial(link="logit")) and 1b.
>
lmer(Y~X+(1|Subject),family=binomial(link="logit"),method="PQL")
>
> Why does 1b produce results different from 1a? The reason why I am
> asking is that the help states that "PQL" is the default of
GLMMs
>
> and
>
> 2.
gamm(Y~X,family=binomial(link="logit"),random=list(Subject=~1))
>
> The interesting thing about the example below is, that gamm is also
> supposed to fit by "PQL". Interestingly, however, the GAMM fit
yields
> about the coefficient estimates of 1b. But the significance values of
> 1a. Any insight would be greatly appreciated.
>
>
> library(lme4)
> library(mgcv)
>
> Y=c(0,1,1,1,1,0,0,0,0,0,1,1,1,1,0,0,0,1,1,1,1)
> X=c(1,2,3,4,3,1,0,0,2,3,3,2,4,3,2,1,1,3,4,2,3)
> Subject=as.factor(c(1,2,3,4,5,6,7,1,2,3,4,5,6,7,1,2,3,4,5,6,7))
> cbind(Y,X,Subject)
>
> r1=lmer(Y~X+(1|Subject),family=binomial(link="logit"))
> summary(r1)
>
>
r2=lmer(Y~X+(1|Subject),family=binomial(link="logit"),method="PQL")
> summary(r2)
>
>
r3=gamm(Y~X,family=binomial(link="logit"),random=list(Subject=~1))
> summary(r3$gam)
>
>
>
> -------------------------
> cuncta stricte discussurus
>
> ______________________________________________
> 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.
>