Frank E Harrell Jr
2004-May-27 06:44 UTC
[R] Getting the same values of adjusted mean and standard errors as SAS
On Thu, 27 May 2004 16:34:58 +0930 "David J. Netherway" <david.netherway at adelaide.edu.au> wrote:> Hello, > > I am trying to get the same values for the adjusted means and standard > errors using R that are given in SAS for the > following data. The model is Measurement ~ Age + Gender + Group. I can > get the adusted means at the mean age > by using predict. I do not know how to get the appropriate standard > errors at the adjusted means for Gender > using values from predict. So I attempted to get them directly from the > residuals as follows. The data is at the end > of the email. While there is a match for the males there is a large > difference for the females indicating that what I am doing is wrong. > > # > meanAge <- mean(dd$Age) > meanAgeM <- mean(dd$Age[d$Gender=="M"]) > meanAgeF <- mean(dd$Age[d$Gender=="F"]). . . . By using sex-specific means of age you are not getting adjusted estimates in the usual sense. I prefer to think of effects as differences in predicted values rather than as complex SAS-like contrasts. The Design package's contrast function makes this easy (including SEs and confidence limits): library(Design) # also requires Hmisc d <- datadist(dd); options(datadist='d') f <- ols(y ~ age + sex + group, data=dd) contrast(f, list(sex='M'), list(sex='F')) # usual adjusted difference M vs F contrast(f, list(sex='M',age=mean(dd$age[dd$sex=='M']), list(sex='F',age=mean(dd$age[dd$sex=='F')) # M vs F not holding age constant You can also experiment with specifying age=tapply(age, sex, mean, na.rm=TRUE) using some of the contrast.Design options. --- Frank E Harrell Jr Professor and Chair School of Medicine Department of Biostatistics Vanderbilt University
David J. Netherway
2004-May-27 07:04 UTC
[R] Getting the same values of adjusted mean and standard errors as SAS
Hello,
I am trying to get the same values for the adjusted means and standard
errors using R that are given in SAS for the
following data. The model is Measurement ~ Age + Gender + Group. I can
get the adusted means at the mean age
by using predict. I do not know how to get the appropriate standard
errors at the adjusted means for Gender
using values from predict. So I attempted to get them directly from the
residuals as follows. The data is at the end
of the email. While there is a match for the males there is a large
difference for the females indicating that what I am doing is wrong.
#
meanAge <- mean(dd$Age)
meanAgeM <- mean(dd$Age[d$Gender=="M"])
meanAgeF <- mean(dd$Age[d$Gender=="F"])
# determine adjusted means for the males and females at meanAge using
predict
# set up data frame to get predicted values at meanAge
evalDF <- data.frame(Age = meanAge, Gender = c("F","M"),
Group = c("1NC","1NC", "2UCLP",
"2UCLP", "3BCLP", "3BCLP",
"4ICP", "4ICP", "5CLPP", "5CLPP"))
mod <-lm(Measurement ~ Age + Gender + Group, data=dd)
pred <- predict(mod, evalDF, se.fit = TRUE)
adjDF <- data.frame(evalDF,fit = pred$fit, se = pred$se.fit)
adjMeanMale <- mean(adjDF$fit[adjDF$Gender=="M"]); # match:
3.889965 cf
3.88996483
adjMeanFemale <- mean(adjDF$fit[adjDF$Gender=="F"]); # match:
3.91111 cf
3.91111036
# Try to get standard errors at the adjusted means for Gender as follows:
ddr <- data.frame(dd, res = residuals(mod))
nM <- summary(ddr$Gender)[["M"]]
seRegM <- sqrt(mean( ddr$res[ddr$Gender=="M"]**2 ))
sxxM <- sum((dd$Age[d$Gender=="M"]-meanAgeM)**2)
syM <- seRegM * sqrt(1/nM + (meanAge-meanAgeM)**2/sxxM); #0.1103335 cf
0.11032602 - matches to 5 decimal places
nF <- summary(ddr$Gender)[["F"]]
seRegF <- sqrt(mean( ddr$res[ddr$Gender=="F"]**2 ))
sxxF <- sum((dd$Age[d$Gender=="F"]-meanAgeF)**2)
syF <- seRegF * sqrt(1/nF + (meanAge-meanAgeF)**2/sxxF); # wrong:
0.07279221 cf 0.14256466
> dd
Measurement Age Gender Group
1 3.8 94 M 3BCLP
2 2.7 88 F 3BCLP
3 3.0 155 M 3BCLP
4 2.7 33 M 3BCLP
5 4.6 109 M 5CLPP
6 5.1 325 M 5CLPP
7 3.9 79 M 5CLPP
8 4.2 126 M 5CLPP
9 3.9 77 F 5CLPP
10 4.0 61 F 5CLPP
11 3.6 49 F 5CLPP
12 3.7 14 F 4ICP
13 4.2 160 F 4ICP
14 3.9 60 M 4ICP
15 5.0 61 M 4ICP
16 3.9 222 F 4ICP
17 3.8 82 F 4ICP
18 4.8 340 F 4ICP
19 3.2 206 M 4ICP
20 3.8 19 M 1NC
21 4.9 166 M 1NC
22 3.8 93 M 1NC
23 3.6 142 M 1NC
24 4.8 241 M 1NC
25 3.9 81 M 1NC
26 4.5 41 M 1NC
27 5.1 244 F 1NC
28 4.6 100 M 1NC
29 5.1 122 F 1NC
30 4.7 194 F 1NC
31 5.1 297 M 1NC
32 3.9 69 M 2UCLP
33 2.5 141 M 2UCLP
34 3.2 104 M 2UCLP
35 3.8 90 M 2UCLP
36 3.8 92 M 2UCLP
37 3.6 149 F 2UCLP
38 3.8 53 F 2UCLP
39 4.7 111 M 2UCLP
40 3.8 116 F 2UCLP
41 3.3 81 M 2UCLP
>