[R] Standard errors from glm

Peter Alspach PAlspach at hortresearch.co.nz
Tue Aug 3 03:10:02 CEST 2004


Roger

Thanks - I have tried that but it doesn't give the standard errors I required.

Using my example:

coef(summary(temp.lm)) gives:

                  Estimate Std. Error     t value     Pr(>|t|)
(Intercept)           44.5  10.535912  4.22364968 0.0003224616
rep2                  -1.0   4.214365 -0.23728369 0.8145375583
trt11                -13.0  14.598987 -0.89047272 0.3824325293
trt12                  3.0  14.598987  0.20549370 0.8389943428
trt13                -17.0  14.598987 -1.16446432 0.2561726432
.. etc ...

that is, standard error for the (intercept) is 10.54, rep 4.21, main effects 14.60, 2-way interactions 20.65 and 3-way interaction 29.20.  These can also be obtained via
sqrt(diag(vcov(temp.lm))).

It is not clear to me how to go from these estimates to those from the aov() call.  I have tried pre- and post- multiplying vcov() by the design matrix but this gives the same standard errors as predict(temp.lm, se=T); i.e. those of the predicted values.

Regards ........

Peter Alspach


>>> "Roger D. Peng" <rpeng at jhsph.edu> 03/08/04 12:22:44 >>>
Try summary(glm.object)$coefficients.

-roger

Peter Alspach wrote:
> Kia ora list members:
> 
> 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.  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.
> 
> predict(lm(S~rep+trt1*trt2*trt3, data=dummy.data),
> type='terms', se=T) or predict(glm(cbind(S,
> 100-S)~rep+trt1*trt2*trt3, data=dummy.data,
> family='binomial'), type='terms', se=T) or, as in my case, 
> predict(glm(cbind(S, 100-S)~rep+trt1*trt2*trt3,
> data=dummy.data, family='quasibinomial'), type='terms', se=T)
> 
> 
> Thanks ........
> 
> Peter Alspach HortResearch
> 
> dummy.data trt1	trt2	trt3	rep	S 0	0	0	1	33 0	0	0	2	55 0	0	1	1
> 18 0	0	1	2	12 0	1	0	1	47 0	1	0	2	16 0	1	1	1	22 0	1	1	2	33 0	2
> 0	1	22 0	2	0	2	18 0	2	1	1	60 0	2	1	2	40 1	0	0	1	38 1	0	0	2	24 
> 1	0	1	1	 8 1	0	1	2	14 1	1	0	1	69 1	1	0	2	42 1	1	1	1	42 1	1	1	2
> 44 1	2	0	1	48 1	2	0	2	26 1	2	1	1	46 1	2	1	2	33 2	0	0	1	48 2	0
> 0	2	46 2	0	1	1	24 2	0	1	2	 8 2	1	0	1	69 2	1	0	2	33 2	1	1	1	22 
> 2	1	1	2	33 2	2	0	1	33 2	2	0	2	18 2	2	1	1	26 2	2	1	2	42 3	0	0	1
> 12 3	0	0	2	42 3	0	1	1	16 3	0	1	2	22 3	1	0	1	14 3	1	0	2	60 3	1
> 1	1	40 3	1	1	2	55 3	2	0	1	47 3	2	0	2	38 3	2	1	1	18 3	2	1	2	44
> 
> 
> ______________________________________________________
> 
> The contents of this e-mail are privileged and/or
> confidenti...{{dropped}}
> 
> ______________________________________________ 
> R-help at stat.math.ethz.ch mailing list 
> https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE
> do read the posting guide!
> http://www.R-project.org/posting-guide.html 
> 

______________________________________________________

The contents of this e-mail are privileged and/or confidenti...{{dropped}}




More information about the R-help mailing list