[R] Standard errors from glm
Prof Brian Ripley
ripley at stats.ox.ac.uk
Tue Aug 3 08:10:56 CEST 2004
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
More information about the R-help
mailing list