> On 04 Aug 2016, at 22:00 , Laviolette, Michael > <michael.laviole...@dhhs.nh.gov> wrote: > > Thanks. I came to the same realization after the original post and was able > to get the correct results with the coefficient vector and covariance matrix > by setting up as a contrast. Linear model contrast estimation in R doesn't > seem straightforward. I turned up several packages, but any recommendations > you might have would be very useful. Thanks again. --M.L.
For full generality contrasts, it doesn't really get much simpler that what you do below. Any "smart" way of specifying your "d" vector ends up with some "d" not specifiable. The only neat thing I can think of is that you can do multiple contrasts with a matrix D instead of a vector D as z <- D %*% coef(fit) and then, although correct, you do not want to get se.fit as sqrt(diag(D %*% vcov(fit) %*% t(D))) if D has a substantial number of rows. Rather, you do it as rowSums((D %*% vcov(fit)) * D). -ps > > a <- 25 # age for which to compute OR's > d <- c(1, 1, a, a) - c(1, 0, a, 0) # contrast > est.ln.or <- crossprod(coef(fit3.16), d) > se.ln.or <- sqrt(t(d) %*% vcov(fit3.16) %*% d) > exp(est.ln.or) > # [1,] 3.899427 > exp(est.ln.or - 1.96 * se.ln.or) > # [1,] 1.712885 > exp(est.ln.or + 1.96 * se.ln.or) > # [1,] 8.877148 > > -----Original Message----- > From: peter dalgaard [mailto:pda...@gmail.com] > Sent: Thursday, August 04, 2016 9:12 AM > To: Laviolette, Michael > Cc: r-help@r-project.org > Subject: Re: [R] Odds ratios in logistic regression models with interaction > > I suspect that "you can't get there from here"... a1$fit and a2$fit are not > independent, so you can't work out the s.e. of their difference using > sqrt(a1$se.fit^2+a2$se.fit^2). > > You need to backtrack a bit and figure out how a1$fit-a2$fit relates to > coef(fit3.16). I suspect it is actually just the age times the interaction > term, but since you give no output and your code uses a bunch of stuff that I > haven't got installed, I can't be bothered to check.... > > Once you have your desired value in the form t(a) %*% coef(...), then use the > result that V(t(a) %*% betahat) == t(a) %*% vcov() %*% a (asymptotically). > > -pd > > On 03 Aug 2016, at 15:08 , Laviolette, Michael > <michael.laviole...@dhhs.nh.gov> wrote: > >> I'm trying to reproduce some results from Hosmer & Lemeshow's "Applied >> Logistic Regression" second edition, pp. 74-79. The objective is to estimate >> odds ratios for low weight births with interaction between mother's age and >> weight (dichotomized at 110 lb.). I can get the point estimates, but I can't >> find an interval option. Can anyone provide guidance on computing the >> confidence intervals? Thanks. -Mike L. >> >> library(dplyr) >> data(birthwt, package = "MASS") >> birthwt <- birthwt %>% >> mutate(low = factor(low, 0:1, c("Not low", "Low")), >> lwd = cut(lwt, c(0, 110, Inf), right = FALSE, >> labels = c("Less than 110", "At least 110")), >> lwd = relevel(lwd, "At least 110")) >> >> # p. 77, Table 3.16, Model 3 >> fit3.16 <- glm(low ~ lwd * age, binomial, birthwt) # p. 78, >> interaction plot visreg::visreg(fit3.16, "age", by = "lwd", xlab = >> "Age", >> ylab = "Estimated logit") # p. 78, covariance matrix >> vcov(fit3.16) >> # odds ratios for ages 15, 20, 25, 30 >> age0 <- seq(15, 30, 5) >> df1 <- data.frame(lwd = "Less than 110", age = age0) >> df2 <- data.frame(lwd = "At least 110", age = age0) >> a1 <- predict(fit3.16, df1, se.fit = TRUE) >> a2 <- predict(fit3.16, df2, se.fit = TRUE) # p. 79, point estimates >> exp(a1$fit - a2$fit) >> >> # How to get CI's? >> # Age OR (95% CI) >> # ---------------------- >> # 15 1.04 (0.29, 3.79) >> # 20 2.01 (0.91, 4.44) >> # 25 3.90 (1.71, 8.88) >> # 30 7.55 (1.95, 29.19) >> >> >> >> >> [[alternative HTML version deleted]] >> >> ______________________________________________ >> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see >> 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. > > -- > Peter Dalgaard, Professor, > Center for Statistics, Copenhagen Business School Solbjerg Plads 3, 2000 > Frederiksberg, Denmark > Phone: (+45)38153501 > Office: A 4.23 > Email: pd....@cbs.dk Priv: pda...@gmail.com > -- Peter Dalgaard, Professor, Center for Statistics, Copenhagen Business School Solbjerg Plads 3, 2000 Frederiksberg, Denmark Phone: (+45)38153501 Office: A 4.23 Email: pd....@cbs.dk Priv: pda...@gmail.com ______________________________________________ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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.
