Solved – Calculating a risk ratio for specific x values from a GAM model using the mgcv package

confidence intervalgeneralized-additive-modelmgcvrratio

I have fit a generalized additive model (GAM) using the mgcv package in R. My model has a dichotomous response variable and so i've used the binomial family link function. After creating the model I would like to do a little post-estimation inference above and beyond the plot.gam graphs.

I would like to take two x-values, for example, and calculate the risk ratio and 95% confidence intervals for that ratio. Obtaining the risk ratio seems fairly straightforward. I could transform the predictions into probabilities and simply divide the two probabilities corresponding to the x-values of interest in order to get the risk ratio. I am less certain how to get the confidence intervals.

In this link here: http://grokbase.com/t/r/r-help/125qbnw21a/r-mgcv-how-to-calculate-a-confidence-interval-of-a-ratio Simon Wood, the author of the mgcv package explained how to get the CIs for a log ratio using a poisson model. I'm uncertain how I would need to change the code to get the risk ratios and 95% CIs from my logistic model.

Here is a reproducible example provided by Simon Wood in the link above:

    library(mgcv)

    ## simulate some data
    dat <- gamSim(1, n=1000, dist="poisson", scale=.25)

    ## fit log-linear model...
    b <- gam(y~s(x0)+s(x1)+s(x2)+s(x3), family=poisson,
    data=dat, method="REML")

    ## data at which predictions to be compared...
    pd <- data.frame(x0=c(.2,.3),x1=c(.5,.5),x2=c(.5,.5),
    x3=c(.5,.5))

    ## log(E(y_1)/E(y_2)) = s(x_1) - s(x_2)
    Xp <- predict(b,newdata=pd,type="lpmatrix")

    ## ... Xp%*%coef(b) gives log(E(y_1)) and log(E(y_2)),
    ## so the required difference is computed as...
    diff <- (Xp[1,]-Xp[2,])
    dly <- t(diff)%*%coef(b) ## required log ratio (diff of logs)
    se.dly <- sqrt(t(diff)%*%vcov(b)%*%diff) ## corresponding s.e.
    dly + c(-2,2)*se.dly ## 95%CI

Any help is greatly appreciated.

Best Answer

This doesn't exactly answer your question, but it might still solve your problem of needing to calculate risk ratios. The epiR package allows you to calculate risk ratios.

I could not get your example to work (see my comment to your question), so here is an example from the package's documentation:

library(epiR) # Used for Risk ratio
library(MASS) # Used for data

dat1 <- birthwt; head(dat1)

## Generate a table of cell frequencies. First set the levels of the outcome
## and the exposure so the frequencies in the 2 by 2 table come out in the
## conventional format:
dat1$low <- factor(dat1$low, levels = c(1,0))
dat1$smoke <- factor(dat1$smoke, levels = c(1,0))
dat1$race <- factor(dat1$race, levels = c(1,2,3))
## Generate the 2 by 2 table. Exposure (rows) = smoke. Outcome (columns) = low.
tab1 <- table(dat1$smoke, dat1$low, dnn = c("Smoke", "Low BW"))
print(tab1)
## Compute the incidence risk ratio and other measures of association:
epi.2by2(dat = tab1, method = "cohort.count", 
conf.level = 0.95, units = 100, outcome = "as.columns")
Related Question