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}}
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 >
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}}
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
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()$contrastsunordered 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}}