[R] Standard errors from glm

Peter Alspach PAlspach at hortresearch.co.nz
Tue Aug 3 01:51:14 CEST 2004


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}}




More information about the R-help mailing list