On Sat, Jan 29, 2011 at 1:05 PM, David Winsemius <dwinsem...@comcast.net> wrote: > > On Jan 29, 2011, at 3:49 PM, Joshua Wiley wrote: > >> On Sat, Jan 29, 2011 at 12:31 PM, David Winsemius >> <dwinsem...@comcast.net> wrote: >>> >>> Huh?. With the same model and data, they should be the same: >> >> I must be missing something, because that is what I would have >> expected too, but it is not what I get (at least when I run the code >> as shown below). >> >>>> lm.mod2 <- lm(out ~ scale(xxA)*scale(xxB), data=dat) >>>> newdata2 <- expand.grid(xxA=c(-1,0,1),xxB=c(-1,0,1)) >>>> newdata2$preds <- predict(lm.mod, newdata2) >> >> ^^^^^ is that lm.mod2? > > UnnHHH . lm.mod2 was the model from Zahn's posting. > > lm.mod <- lm(out ~ I(scale(xxA))*I(scale(xxB)), data=dat) > > So I guess you need the "I" to get the correct results with the scaled > numbers, but this "works" with the unscaled numbers, although the numbers > (meanA, meanB) are somewhat different, which deserves investigation. The > pair (meanA, meanB) gives the same result but the others don't seem to get > to the same numbers. Which is correct? > >> meanA<-mean(dat$xxA); meanB<-mean(dat$xxB); sdA<-sd(dat$xxA); >> sdB<-sd(dat$xxB) >> newdata2 <- >> expand.grid(xxA=c(meanA-sdA,meanA,meanA+sdA),xxB=c(meanB-sdB,meanB,meanB+sdB)) >> newdata2$preds <- predict(lm.mod2, newdata2) >> newdata2 > xxA xxB preds > 1 9.304820 18.75866 222.1326 > 2 9.905043 18.75866 234.0992 > 3 10.505266 18.75866 246.0659 > 4 9.304820 19.73195 232.7828 > 5 9.905043 19.73195 244.8696 > 6 10.505266 19.73195 256.9565 > 7 9.304820 20.70524 243.4330 > 8 9.905043 20.70524 255.6400 > 9 10.505266 20.70524 267.8471 >
I may be completely barking up the wrong tree, but could it have something to do with the interaction term being created from scaled vs. unscaled terms? > xxAB <- dat$xxA * dat$xxB > meanAB <- meanA * meanB > sdAB <- sqrt(sdA^2 + sdB^2 + meanA^2 * sdB^2 + meanB^2 * sdA^2) > (xxAB - meanAB)/sdAB [1] 0.17672311 -0.93203381 -1.39257245 -1.31620813 0.08964000 0.71663046 [7] -1.21984785 -0.39528381 -1.40027554 0.06304763 0.57713584 1.30481682 [13] 0.76070783 1.73067640 0.49481401 -0.22607119 -0.40356164 -0.70588375 [19] 1.34382743 0.65508728 > (sxxAB <- with(dat, scale(xxA) * scale(xxB))[1:20,]) [1] 0.01838999 0.25899229 -0.15337178 1.07536361 -0.26520706 0.25119036 [7] -0.11295999 0.05205466 -1.49676116 -0.14367368 -1.86921225 0.80860167 [13] -0.44816896 1.44122535 -0.73361472 -0.14930187 -1.46192081 0.19784309 [19] 0.62089372 0.0821291515 > >>>> newdata2 >>> >>> xxA xxB preds >>> 1 -1 -1 218.6366 >>> 2 0 -1 232.4330 >>> 3 1 -1 246.2295 >>> 4 -1 0 230.9129 >>> 5 0 0 244.8696 >>> 6 1 0 258.8263 >>> 7 -1 1 243.1892 >>> 8 0 1 257.3062 >>> 9 1 1 271.4232 >>> >>>> Predict(rms.res,xxA=c(-1,0,1),xxB=c(-1,0,1), conf.int=FALSE) >>> >>> xxA xxB yhat >>> 1 -1 -1 218.6366 >>> 2 0 -1 232.4330 >>> 3 1 -1 246.2295 >>> 4 -1 0 230.9129 >>> 5 0 0 244.8696 >>> 6 1 0 258.8263 >>> 7 -1 1 243.1892 >>> 8 0 1 257.3062 >>> 9 1 1 271.4232 >> >> This is copied directly from my console in a clean session. I get >> different values for the predicted values from predict(glm or lm) and >> "yhat" from Predict(Glm). >> >> >>> sessionInfo() >> >> R version 2.12.1 (2010-12-16) >> Platform: i486-pc-linux-gnu (32-bit) >> >> locale: >> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C >> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 >> [5] LC_MONETARY=C LC_MESSAGES=en_US.UTF-8 >> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C >> [9] LC_ADDRESS=C LC_TELEPHONE=C >> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C >> >> attached base packages: >> [1] stats graphics grDevices utils datasets methods base >>> >>> ################################# >>> require(rms) >> >> Loading required package: rms >> Loading required package: Hmisc >> Loading required package: survival >> Loading required package: splines >> >> Attaching package: 'Hmisc' >> >> The following object(s) are masked from 'package:survival': >> >> untangle.specials >> >> The following object(s) are masked from 'package:base': >> >> format.pval, round.POSIXt, trunc.POSIXt, units >> >> >> Attaching package: 'rms' >> >> The following object(s) are masked from 'package:survival': >> >> Surv >> >>> set.seed(10) >>> dat <- data.frame(xxA = rnorm(20, 10), xxB = rnorm(20, 20)) >>> dat$out <- with(dat, xxA+xxB+xxA*xxB+rnorm(20,20)) >>> >>> rms.res <- Glm(out ~ scale(xxA) * scale(xxB), data=dat) >>> glm.mod <- glm(out ~ scale(xxA) * scale(xxB), data=dat) >>> lm.mod2 <- lm(out ~ scale(xxA)*scale(xxB), data=dat) >>> >>> newdata2 <- expand.grid(xxA = -1:1, xxB = -1:1) >>> newdata2$preds1 <- predict(glm.mod, newdata2) >>> newdata2$preds2 <- predict(lm.mod2, newdata2) >>> newdata2 >> >> xxA xxB preds1 preds2 >> 1 -1 -1 143.3433 143.3433 >> 2 0 -1 131.6110 131.6110 >> 3 1 -1 119.8787 119.8787 >> 4 -1 0 137.1149 137.1149 >> 5 0 0 126.9708 126.9708 >> 6 1 0 116.8267 116.8267 >> 7 -1 1 130.8865 130.8865 >> 8 0 1 122.3306 122.3306 >> 9 1 1 113.7747 113.7747 >>> >>> Predict(rms.res, xxA= -1:1, xxB= -1:1, conf.int=FALSE) >> >> xxA xxB yhat >> 1 -1 -1 212.2324 >> 2 0 -1 229.6472 >> 3 1 -1 247.0620 >> 4 -1 0 221.7795 >> 5 0 0 240.6413 >> 6 1 0 259.5030 >> 7 -1 1 231.3266 >> 8 0 1 251.6353 >> 9 1 1 271.9441 >> >> Response variable (y): X * Beta >>> >>> ################################## > > David Winsemius, MD > West Hartford, CT > > -- Joshua Wiley Ph.D. Student, Health Psychology University of California, Los Angeles http://www.joshuawiley.com/ ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
