[R] Odds ratios in logistic regression models with interaction

[peter dalgaard]

> On 04 Aug 2016, at 22:00 , Laviolette, Michael <Michael.Laviolette at 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: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
> 

-- 
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