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" <[EMAIL PROTECTED]> 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}}
>
> ______________________________________________
> [EMAIL PROTECTED] 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}}
______________________________________________
[EMAIL PROTECTED] 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
