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 >