R help - Jan 2006 - multiple GLMs with lmList in lme4

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

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
>
> ______________________________________________
> R-help at stat.math.ethz.ch 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
tlumley at u.washington.edu	University of Washington, Seattle

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 <tlumley at u.washington.edu> 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
>
> ______________________________________________
> R-help at stat.math.ethz.ch 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
tlumley at u.washington.edu	University of Washington, Seattle