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

Question

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.

Answer 1

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