[R] Standard errors from glm
Peter Alspach
PAlspach at hortresearch.co.nz
Wed Aug 4 02:23:15 CEST 2004
Dear Prof Ripley
Thanks for your reply and clarification. However:
1. Regarding model.tables() returning "Design is unbalanced". Setting contrasts to Helmert does indeed make the design balanced, but model.tables() still returns "Design is unbalanced":
> options()$contrasts
unordered ordered
"contr.treatment" "contr.poly"
> aov(S~rep+trt1*trt2*trt3, data=dummy.data)
Call:
...
Residual standard error: 14.59899
Estimated effects may be unbalanced
> options(contrasts=c("contr.helmert", "contr.treatment"))
> aov(S~rep+trt1*trt2*trt3, data=dummy.data)
Call:
...
Residual standard error: 14.59899
Estimated effects are balanced
> model.tables(aov(S~rep+trt1*trt2*trt3, data=dummy.data), se=T)
Design is unbalanced - use se.contrasts for se's
Tables of effects
...
However, this is a relatively minor issue, and covered in ?model.tables which clearly states that "The implementation is incomplete, and only the simpler cases have been tested thoroughly."
2. You point out that "In either case you can predict something you want to estimate and use
predict(, se=TRUE)." Doesn't this give the standard error of the predicted value, rather than the mean for, say, trt1 level 0? For example:
> predict(temp.lm, newdata=data.frame(rep='1', trt1='0', trt2='1', trt3='0'), se=T)
$fit
[1] 32
$se.fit
[1] 10.53591
$df
[1] 23
$residual.scale
[1] 14.59899
Whereas from the analysis of variance table we can get the standard error of the mean for trt1 as being sqrt(anova(temp.lm)[9,3]/12) = 4.214365. It is the equivalent of this latter value that I'm after in the glm() case.
>>> Prof Brian Ripley <ripley at stats.ox.ac.uk> 03/08/04 18:10:56 >>>
On Tue, 3 Aug 2004, Peter Alspach wrote:
[Lines wrapped for legibility.]
> I'm having a little difficulty getting the correct standard errors from
> a glm.object (R 1.9.0 under Windows XP 5.1). predict() will gives
> standard errors of the predicted values, but I am wanting the standard
> errors of the mean.
>
> To clarify:
>
> Assume I have a 4x3x2 factorial with 2 complete replications (i.e. 48
> observations, I've appended a dummy set of data at the end of this
> message). Call the treatments trt1 (4 levels), trt2 (3 levels) and trt3
> (2 levels) and the replications rep - all are factors. The observed
> data is S. Then:
>
> temp.aov <- aov(S~rep+trt1*trt2*trt3, data=dummy.data)
> model.tables(temp.aov, type='mean', se=T)
>
> Returns the means, but states "Design is unbalanced - use se.contrasts
> for se's" which is a little surprising since the design is balanced.
If you used the default treatment contrasts, it is not. Try Helmert
contrasts with aov().
> Nevertheless, se.contrast gives what I'd expect:
>
> se.contrast(temp.aov, list(trt1==0, trt1==1), data=dummy.data)
> [1] 5.960012
>
> i.e. standard error of mean is 5.960012/sqrt(2) = 4.214, which is the
> sqrt(anova(temp.aov)[9,3]/12) as expected. Similarly for interactions,
> e.g.:
>
> se.contrast(temp.aov, list(trt1==0 & trt2==0, trt1==1 & trt2==1), data=dummy.data)/sqrt(2)
> [1] 7.299494
>
> How do I get the equivalent of these standard errors if I have used
> lm(), and by extension glm()? I think I should be able to get these
> using predict(..., type='terms', se=T) or coef(summary()) but can't
> quite see how.
In either case you can predict something you want to estimate and use
predict(, se=TRUE).
--
Brian D. Ripley, ripley at stats.ox.ac.uk
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UK Fax: +44 1865 272595
______________________________________________________
The contents of this e-mail are privileged and/or confidenti...{{dropped}}
More information about the R-help
mailing list