Odds ratios for interaction between continuous variables

[jCeradini]

I'm trying to work through the formula for the 95% CI for the OR but the estimates I'm getting don't make sense (given my probability plots at the top of this thread and the fact that the interaction is marginally significant). Am I misinterpreting the variance formula? I've tripled checked my code and haven't found any syntax errors....

Thanks again.
Joe

Output from Huggins robust design:

Code: Select all

estimate        se        lcl        ucl          var           pred
-0.7805226 0.7506744 -2.2518445  0.6907993 0.5635120548  S:(Intercept)
-4.6189001 2.0282830 -8.5943350 -0.6434653 4.1139319281  S:grass
0.0120384  0.0295291 -0.0458387  0.0699155 0.0008719677  S:shrub
0.7729132  0.5820680 -0.3679401  1.9137665 0.3388031566  S:season2014
0.1835990  0.0894282  0.0083198  0.3588782 0.0079974030  S:grass:shrub

Calculate odds ratio and 95% CI for interaction, grass*shrub:

Code: Select all

# Different values of shrub cover
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
s.cov0 <- 0; s.cov15 <- 15; s.cov30 <- 30

# 1) Mean log(odds) for grass effect at constant shrub value
# Formula: grass beta + (interaction beta * shrub cover value)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
grass.shrub.0 <- -4.6189001 + (0.1835990*s.cov0)
grass.shrub.15 <- -4.6189001 + (0.1835990*s.cov15)
grass.shrub.30 <- -4.6189001 + (0.1835990*s.cov30)

# 2) Variance of betas:
#~~~~~~~~~~~~~~~~~~~~~~
top.betas <- model$results$beta
top.betas$var <- top.betas$se^2
top.betas$pred <- row.names(top.betas)
grass.var <- top.betas$var[top.betas$pred == "grass"]
grassXshrub.var <- top.betas$var[top.betas$pred == "grass:shrub"]

# 3) Covariance between grass beta (2) and interaction beta (5)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
grass.grassXshrub.cov <- model$results$beta.vcv[2, 5]

# 4) Calculate 95% CI for OR at different shrub values
# Shrub = 15
#~~~~~~~~~~~
grass.shrub.15.var <- grass.var + ((s.cov15^2)*grassXshrub.var) - ((2*s.cov15)*grass.grassXshrub.cov)

log.CI.15 <- 1.96 * sqrt(grass.shrub.15.var)
exp(grass.shrub.15 - log.CI.15) # 0.000227
exp(grass.shrub.15 + log.CI.15) # 105.501

# Shrub = 30
#~~~~~~~~~~~~
grass.shrub.30.var <- grass.var + ((s.cov30^2)*grassXshrub.var) - ((2*s.cov30)*grass.grassXshrub.cov)

log.CI.30 <- 1.96 * sqrt(grass.shrub.30.var)
exp(grass.shrub.30 - log.CI.30) # 0.00026
exp(grass.shrub.30 + log.CI.30) # 22181.41

[jlaake]

Sorry. I just looked at my post which had a typo. It should have been +2*s*cov() rather than -.

[jCeradini]

No problem - thanks for double checking.

[cooch]

Sorry. I just looked at my post which had a typo. It should have been +2*s*cov() rather than -.

Missed that too, although it is correct in the book -- example on pp. 78-80, loosely drawn from an example similar to the original post in this thread. Obviously, getting the sign right can make a big difference if sampling covariance is > negligible.

[jCeradini]

Example of why you should check for updates of the MARK book often:
I read the effect size sections in ch 6 of the MARK book before posting this question but I'm using an outdated copy of the MARK book. My section on beta estimates and odds ratios stops after the good/poor colony example. Thanks to Jeff and Evan, I now see that the section actually continues in the newer version with an example on calculating odds ratios with continuous variables (p. 78 - 80) that would have saved us all a lot of emails....my mistake.

Maybe some useful R code with a loop for estimating OR and CI's at different covariate values when dealing with an interaction:

Code: Select all

# s = shrub covariate value, B_1 = grass beta, B_3 = interaction beta

# FORMULA 1: MEAN log(Odds) for one unit change
# B_1 + B_3s

# FORMULA 2: variance for the mean log(odds)
# var(log(O))= var(B_1) + s^2var(B_3) + 2s*cov(B_1,B_3)

# FORMULA 3: 95% CI for mean OR and backtransformation to OR scale
# LCL.log = log(O) - 1.96*sqrt{var(log(O)}
# UCL.log = log(O) + 1.96*sqrt{var(log(O)}
# Back-transform CI's to OR scale:
# exp(LCL.log), exp(UCL.log)

# Variance of betas:
#~~~~~~~~~~~~~~~~~~~~~
top.betas <- model$results$beta
top.betas$var <- top.betas$se^2
top.betas$pred <- row.names(top.betas)

grass.var <- top.betas$var[top.betas$pred == "S:grass"]
grassXshrub.var <- top.betas$var[top.betas$pred == "S:grass:shrub"]

# Covariance between grass beta (2) and interaction (5)
grass.grassXshrub.cov <- model$results$beta.vcv[2, 5]


# Loop
#~~~~~~~~
shrub.val <- seq(0, 30, 2)
grass.OR <- data.frame(shrub.val = seq(0, 30, 2), mean.OR = rep(NA, length(shrub.val)), lcl.OR = rep(NA, length(shrub.val)), ucl.OR = rep(NA, length(shrub.val)))

for(i in 1:length(shrub.val)){
  # Mean log(odds) for grass
  temp.grass.mean <- -4.6189001 + (0.1835990*shrub.val[i])
 
  # Variance of log(odds) for grass effect
  temp.var <- grass.var + ((shrub.val[i]^2)*grassXshrub.var) +
                      ((2*shrub.val[i])*grass.grassXshrub.cov)

  # Log 95% CI
  temp.CI <- 1.96 * sqrt(temp.var)
 
  # Mean and 95% CI for odds ratio (back-transformed)
  grass.OR$lcl.OR[i] <- exp(temp.grass.mean - temp.CI)
  grass.OR$ucl.OR[i] <- exp(temp.grass.mean + temp.CI)
  grass.OR$mean.OR[i] <- exp(temp.grass.mean)
}