Laviolette, Michael
2016-Aug-03 13:08 UTC
[R] Odds ratios in logistic regression models with interaction
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]]
peter dalgaard
2016-Aug-04 13:12 UTC
[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.Laviolette at 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 at 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.mes at cbs.dk Priv: PDalgd at gmail.com
Laviolette, Michael
2016-Aug-04 20:00 UTC
[R] Odds ratios in logistic regression models with interaction
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. 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:pdalgd at gmail.com] Sent: Thursday, August 04, 2016 9:12 AM To: Laviolette, Michael Cc: r-help at 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.Laviolette at 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 at 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.mes at cbs.dk Priv: PDalgd at gmail.com