extracting coefficients from lmer

classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

extracting coefficients from lmer

Jonathan Williams-2
Dear R-Helpers,

I want to compare the results of outputs from glmmPQL and lmer analyses.
I could do this if I could extract the coefficients and standard errors
from the summaries of the lmer models. This is easy to do for the glmmPQL
summaries, using

> glmm.fit <- try(glmmPQL(score ~ x*type, random = ~ 1 | subject, data = df,
        family = binomial), TRUE)
> summary(glmmPQL.fit)$tTable

Linear mixed-effects model fit by maximum likelihood
 Data: df
       AIC      BIC    logLik
  1800.477 1840.391 -890.2384

Random effects:
 Formula: ~1 | subject
        (Intercept)  Residual
StdDev:   0.6355517 0.9650671

Variance function:
 Structure: fixed weights
 Formula: ~invwt
Fixed effects: score ~ x * type
                 Value Std.Error  DF    t-value p-value
(Intercept) -0.0812834 0.2933314 294 -0.2771043  0.7819
x1           0.4143072 0.4180624  98  0.9910176  0.3241
type2        0.8509166 0.4084443 294  2.0833112  0.0381
type3        0.6691275 0.4024369 294  1.6626894  0.0974
type4       -0.7830413 0.4123851 294 -1.8988109  0.0586
x1:type2     1.0643239 0.6791126 294  1.5672274  0.1181
x1:type3    -0.7533085 0.5674532 294 -1.3275253  0.1854
x1:type4    -0.0549616 0.5777216 294 -0.0951351  0.9243
etc.

However, there seems to be no route to extract the corresponding information
from the lmer model:-

> lmer.fit=try(lmer(score~x*type+(1|subject), data=df, family=binomial,
        method='AGQ'),TRUE)
> summary(lmer.fit)

Generalized linear mixed model fit using AGQ
Formula: score ~ x * type + (1 | subject)
   Data: df
 Family: binomial(logit link)
      AIC      BIC    logLik deviance
 510.2616 550.1762 -245.1308 490.2616
Random effects:
     Groups        Name    Variance    Std.Dev.
    subject (Intercept)     0.46269     0.68021
# of obs: 400, groups: subject, 100

Estimated scale (compare to 1)  1.019134

Fixed effects:
             Estimate Std. Error  z value Pr(>|z|)
(Intercept) -0.087284   0.300896 -0.29008  0.77175
x1           0.446289   0.428844  1.04068  0.29803
type2        0.913571   0.418978  2.18047  0.02922 *
type3        0.719023   0.412816  1.74175  0.08155 .
type4       -0.839842   0.423021 -1.98534  0.04711 *
x1:type2     1.112673   0.696629  1.59722  0.11022
x1:type3    -0.809599   0.582089 -1.39085  0.16427
x1:type4    -0.062235   0.592623 -0.10502  0.91636
etc.

> summary(lmer.fit)$tTable
NULL
> names(summary(lmer.fit))
NULL
> names(lmer.fit)
NULL
> lmer.fit$coef
NULL

So, then I tried to find out if lmer returns different information.
> help(lmer)
This says "see lmer-class" ->
> help(lmer-class)
No documentation for 'lmer - class' in specified packages and libraries:
you could try 'help.search("lmer - class")'

So, then I tried
> help.search('lmer-class')

This returns
Help files with alias or concept or title matching 'lmer-class' using fuzzy
matching:
lmer-class(Matrix)               Mixed model representations

So, I loaded library Matrix and tried again
> library(Matrix)
> help(lmer-class)
No documentation for 'lmer - class' in specified packages and libraries:
you could try 'help.search("lmer - class")'

If someone could tell me how to extract the Estimates and Std. Errors from
the lmer summary, I'd be very grateful. I would also be very grateful if
someone
could let me know why the Estimates and Std. Errors from the lmer model are
both
larger than those from the glmmPQL model.

I am running R 2.2.1 on a Windows XP machine.

Thanks,

Jonathan Williams

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Reply | Threaded
Open this post in threaded view
|

Re: extracting coefficients from lmer

Dimitris Rizopoulos
you can use something like:

Vcov <- vcov(lmer.fit, useScale = FALSE)
betas <- fixef(lmer.fit)
se <- sqrt(diag(Vcov))
zval <- betas / se
pval <- 2 * pnorm(abs(zval), lower.tail = FALSE)
###############
cbind(betas, se, zval, pval)


I hope it helps.

Best,
Dimitris

----
Dimitris Rizopoulos
Ph.D. Student
Biostatistical Centre
School of Public Health
Catholic University of Leuven

Address: Kapucijnenvoer 35, Leuven, Belgium
Tel: +32/(0)16/336899
Fax: +32/(0)16/337015
Web: http://www.med.kuleuven.be/biostat/
     http://www.student.kuleuven.be/~m0390867/dimitris.htm


----- Original Message -----
From: "Jonathan Williams"
<[hidden email]>
To: "Ethz. Ch" <[hidden email]>
Sent: Tuesday, January 10, 2006 11:54 AM
Subject: [R] extracting coefficients from lmer


> Dear R-Helpers,
>
> I want to compare the results of outputs from glmmPQL and lmer
> analyses.
> I could do this if I could extract the coefficients and standard
> errors
> from the summaries of the lmer models. This is easy to do for the
> glmmPQL
> summaries, using
>
>> glmm.fit <- try(glmmPQL(score ~ x*type, random = ~ 1 | subject,
>> data = df,
> family = binomial), TRUE)
>> summary(glmmPQL.fit)$tTable
>
> Linear mixed-effects model fit by maximum likelihood
> Data: df
>       AIC      BIC    logLik
>  1800.477 1840.391 -890.2384
>
> Random effects:
> Formula: ~1 | subject
>        (Intercept)  Residual
> StdDev:   0.6355517 0.9650671
>
> Variance function:
> Structure: fixed weights
> Formula: ~invwt
> Fixed effects: score ~ x * type
>                 Value Std.Error  DF    t-value p-value
> (Intercept) -0.0812834 0.2933314 294 -0.2771043  0.7819
> x1           0.4143072 0.4180624  98  0.9910176  0.3241
> type2        0.8509166 0.4084443 294  2.0833112  0.0381
> type3        0.6691275 0.4024369 294  1.6626894  0.0974
> type4       -0.7830413 0.4123851 294 -1.8988109  0.0586
> x1:type2     1.0643239 0.6791126 294  1.5672274  0.1181
> x1:type3    -0.7533085 0.5674532 294 -1.3275253  0.1854
> x1:type4    -0.0549616 0.5777216 294 -0.0951351  0.9243
> etc.
>
> However, there seems to be no route to extract the corresponding
> information
> from the lmer model:-
>
>> lmer.fit=try(lmer(score~x*type+(1|subject), data=df,
>> family=binomial,
> method='AGQ'),TRUE)
>> summary(lmer.fit)
>
> Generalized linear mixed model fit using AGQ
> Formula: score ~ x * type + (1 | subject)
>   Data: df
> Family: binomial(logit link)
>      AIC      BIC    logLik deviance
> 510.2616 550.1762 -245.1308 490.2616
> Random effects:
>     Groups        Name    Variance    Std.Dev.
>    subject (Intercept)     0.46269     0.68021
> # of obs: 400, groups: subject, 100
>
> Estimated scale (compare to 1)  1.019134
>
> Fixed effects:
>             Estimate Std. Error  z value Pr(>|z|)
> (Intercept) -0.087284   0.300896 -0.29008  0.77175
> x1           0.446289   0.428844  1.04068  0.29803
> type2        0.913571   0.418978  2.18047  0.02922 *
> type3        0.719023   0.412816  1.74175  0.08155 .
> type4       -0.839842   0.423021 -1.98534  0.04711 *
> x1:type2     1.112673   0.696629  1.59722  0.11022
> x1:type3    -0.809599   0.582089 -1.39085  0.16427
> x1:type4    -0.062235   0.592623 -0.10502  0.91636
> etc.
>
>> summary(lmer.fit)$tTable
> NULL
>> names(summary(lmer.fit))
> NULL
>> names(lmer.fit)
> NULL
>> lmer.fit$coef
> NULL
>
> So, then I tried to find out if lmer returns different information.
>> help(lmer)
> This says "see lmer-class" ->
>> help(lmer-class)
> No documentation for 'lmer - class' in specified packages and
> libraries:
> you could try 'help.search("lmer - class")'
>
> So, then I tried
>> help.search('lmer-class')
>
> This returns
> Help files with alias or concept or title matching 'lmer-class'
> using fuzzy
> matching:
> lmer-class(Matrix)               Mixed model representations
>
> So, I loaded library Matrix and tried again
>> library(Matrix)
>> help(lmer-class)
> No documentation for 'lmer - class' in specified packages and
> libraries:
> you could try 'help.search("lmer - class")'
>
> If someone could tell me how to extract the Estimates and Std.
> Errors from
> the lmer summary, I'd be very grateful. I would also be very
> grateful if
> someone
> could let me know why the Estimates and Std. Errors from the lmer
> model are
> both
> larger than those from the glmmPQL model.
>
> I am running R 2.2.1 on a Windows XP machine.
>
> Thanks,
>
> Jonathan Williams
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
> http://www.R-project.org/posting-guide.html
>


Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html