multiple GLMs with lmList in lme4

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

multiple GLMs with lmList in lme4

Daniel Farewell
I'd like to fit a GLM to each of a number of subsets of some data. The `family' argument to `lmList' (in lme4) has given me cause for optimism, but so far I've only been able to achieve linear model fits. For example

> df <- data.frame(gp = gp.temp <- factor(rep(1:3, each = 100)),
x = x.temp <- rnorm(300),
y = rpois(300, exp((-1:1)[gp.temp] + x.temp)))

Then a call to `glm' on the group 1 subset gives

> glm(y ~ x, family = poisson, data = df, subset = gp == 1)

Call:  glm(formula = y ~ x, family = poisson, data = df, subset = gp == 1)

Coefficients:
(Intercept)            x  
    -1.0143       0.9726  

Degrees of Freedom: 99 Total (i.e. Null);  98 Residual
Null Deviance:      138.5
Residual Deviance: 82.76        AIC: 178.5

(the right answer) but `lmList' gives

> show(lmList(y ~ x | gp, family = poisson, data = df))
Call: lmList(formula = y ~ x | gp, data = df, family = poisson)
Coefficients:
  (Intercept)         x
1   0.5560377 0.6362124
2   1.8431794 1.8541193
3   4.5773106 4.7871929

Degrees of freedom: 300 total; 294 residual
Residual standard error: 2.655714

which come from linear model fits, e.g.

> lm(y ~ x, data = df, subset = gp == 1)

Call:
lm(formula = y ~ x, data = df, subset = gp == 1)

Coefficients:
(Intercept)            x  
     0.5560       0.6362

Any suggestions as to why lmList matches the linear fits rather than the GLM fits would be greatly appreciated. I'm using R2.2.1 with lme version 0.98-1 in Windows XP.

Daniel Farewell
Cardiff University

______________________________________________
[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: multiple GLMs with lmList in lme4

Thomas Lumley
On Tue, 17 Jan 2006, Daniel Farewell wrote:

> I'd like to fit a GLM to each of a number of subsets of some data. The
> `family' argument to `lmList' (in lme4) has given me cause for optimism,
> but so far I've only been able to achieve linear model fits. For example
>
>> df <- data.frame(gp = gp.temp <- factor(rep(1:3, each = 100)),
> x = x.temp <- rnorm(300),
> y = rpois(300, exp((-1:1)[gp.temp] + x.temp)))

Unless you are particularly attached to lmList() you might try by():

by(df,df$gp,function(subset) glm(y~x,family=poisson,data=subset))

  -thomas


>
> Then a call to `glm' on the group 1 subset gives
>
>> glm(y ~ x, family = poisson, data = df, subset = gp == 1)
>
> Call:  glm(formula = y ~ x, family = poisson, data = df, subset = gp == 1)
>
> Coefficients:
> (Intercept)            x
>    -1.0143       0.9726
>
> Degrees of Freedom: 99 Total (i.e. Null);  98 Residual
> Null Deviance:      138.5
> Residual Deviance: 82.76        AIC: 178.5
>
> (the right answer) but `lmList' gives
>
>> show(lmList(y ~ x | gp, family = poisson, data = df))
> Call: lmList(formula = y ~ x | gp, data = df, family = poisson)
> Coefficients:
>  (Intercept)         x
> 1   0.5560377 0.6362124
> 2   1.8431794 1.8541193
> 3   4.5773106 4.7871929
>
> Degrees of freedom: 300 total; 294 residual
> Residual standard error: 2.655714
>
> which come from linear model fits, e.g.
>
>> lm(y ~ x, data = df, subset = gp == 1)
>
> Call:
> lm(formula = y ~ x, data = df, subset = gp == 1)
>
> Coefficients:
> (Intercept)            x
>     0.5560       0.6362
>
> Any suggestions as to why lmList matches the linear fits rather than the GLM fits would be greatly appreciated. I'm using R2.2.1 with lme version 0.98-1 in Windows XP.
>
> Daniel Farewell
> Cardiff University
>
> ______________________________________________
> [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
>

Thomas Lumley Assoc. Professor, Biostatistics
[hidden email] University of Washington, Seattle

______________________________________________
[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: multiple GLMs with lmList in lme4

Daniel Farewell
In reply to this post by Daniel Farewell
Many thanks! That'll work great. It's always good to discover a new, general, function like by().

I would still be interested to know how the family argument to lmList() should be used.

Daniel

>>> Thomas Lumley <[hidden email]> 01/17/06 3:15 pm >>>
On Tue, 17 Jan 2006, Daniel Farewell wrote:

> I'd like to fit a GLM to each of a number of subsets of some data. The
> `family' argument to `lmList' (in lme4) has given me cause for optimism,
> but so far I've only been able to achieve linear model fits. For example
>
>> df <- data.frame(gp = gp.temp <- factor(rep(1:3, each = 100)),
> x = x.temp <- rnorm(300),
> y = rpois(300, exp((-1:1)[gp.temp] + x.temp)))

Unless you are particularly attached to lmList() you might try by():

by(df,df$gp,function(subset) glm(y~x,family=poisson,data=subset))

  -thomas


>
> Then a call to `glm' on the group 1 subset gives
>
>> glm(y ~ x, family = poisson, data = df, subset = gp == 1)
>
> Call:  glm(formula = y ~ x, family = poisson, data = df, subset = gp == 1)
>
> Coefficients:
> (Intercept)            x
>    -1.0143       0.9726
>
> Degrees of Freedom: 99 Total (i.e. Null);  98 Residual
> Null Deviance:      138.5
> Residual Deviance: 82.76        AIC: 178.5
>
> (the right answer) but `lmList' gives
>
>> show(lmList(y ~ x | gp, family = poisson, data = df))
> Call: lmList(formula = y ~ x | gp, data = df, family = poisson)
> Coefficients:
>  (Intercept)         x
> 1   0.5560377 0.6362124
> 2   1.8431794 1.8541193
> 3   4.5773106 4.7871929
>
> Degrees of freedom: 300 total; 294 residual
> Residual standard error: 2.655714
>
> which come from linear model fits, e.g.
>
>> lm(y ~ x, data = df, subset = gp == 1)
>
> Call:
> lm(formula = y ~ x, data = df, subset = gp == 1)
>
> Coefficients:
> (Intercept)            x
>     0.5560       0.6362
>
> Any suggestions as to why lmList matches the linear fits rather than the GLM fits would be greatly appreciated. I'm using R2.2.1 with lme version 0.98-1 in Windows XP.
>
> Daniel Farewell
> Cardiff University
>
> ______________________________________________
> [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 
>

Thomas Lumley Assoc. Professor, Biostatistics
[hidden email] University of Washington, Seattle

______________________________________________
[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: multiple GLMs with lmList in lme4

Douglas Bates
On 1/17/06, Daniel Farewell <[hidden email]> wrote:
> Many thanks! That'll work great. It's always good to discover a new, general, function like by().
>
> I would still be interested to know how the family argument to lmList() should be used.

As you have discovered the family argument has no effect in lmList()
from the current version of the lme4 package.   A new version of the
package with this fixed will be uploaded shortly.

> Daniel
>
> >>> Thomas Lumley <[hidden email]> 01/17/06 3:15 pm >>>
> On Tue, 17 Jan 2006, Daniel Farewell wrote:
>
> > I'd like to fit a GLM to each of a number of subsets of some data. The
> > `family' argument to `lmList' (in lme4) has given me cause for optimism,
> > but so far I've only been able to achieve linear model fits. For example
> >
> >> df <- data.frame(gp = gp.temp <- factor(rep(1:3, each = 100)),
> > x = x.temp <- rnorm(300),
> > y = rpois(300, exp((-1:1)[gp.temp] + x.temp)))
>
> Unless you are particularly attached to lmList() you might try by():
>
> by(df,df$gp,function(subset) glm(y~x,family=poisson,data=subset))
>
>         -thomas
>
>
> >
> > Then a call to `glm' on the group 1 subset gives
> >
> >> glm(y ~ x, family = poisson, data = df, subset = gp == 1)
> >
> > Call:  glm(formula = y ~ x, family = poisson, data = df, subset = gp == 1)
> >
> > Coefficients:
> > (Intercept)            x
> >    -1.0143       0.9726
> >
> > Degrees of Freedom: 99 Total (i.e. Null);  98 Residual
> > Null Deviance:      138.5
> > Residual Deviance: 82.76        AIC: 178.5
> >
> > (the right answer) but `lmList' gives
> >
> >> show(lmList(y ~ x | gp, family = poisson, data = df))
> > Call: lmList(formula = y ~ x | gp, data = df, family = poisson)
> > Coefficients:
> >  (Intercept)         x
> > 1   0.5560377 0.6362124
> > 2   1.8431794 1.8541193
> > 3   4.5773106 4.7871929
> >
> > Degrees of freedom: 300 total; 294 residual
> > Residual standard error: 2.655714
> >
> > which come from linear model fits, e.g.
> >
> >> lm(y ~ x, data = df, subset = gp == 1)
> >
> > Call:
> > lm(formula = y ~ x, data = df, subset = gp == 1)
> >
> > Coefficients:
> > (Intercept)            x
> >     0.5560       0.6362
> >
> > Any suggestions as to why lmList matches the linear fits rather than the GLM fits would be greatly appreciated. I'm using R2.2.1 with lme version 0.98-1 in Windows XP.
> >
> > Daniel Farewell
> > Cardiff University
> >
> > ______________________________________________
> > [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
> >
>
> Thomas Lumley                   Assoc. Professor, Biostatistics
> [hidden email]        University of Washington, Seattle
>
> ______________________________________________
> [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
>

______________________________________________
[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