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

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.

## migrated from stackoverflow.comNov 4 '15 at 11:16

This question came from our site for professional and enthusiast programmers.

• Is the @ in gamSim(1,n@0,dist="poisson",scale=.25) a typo? When I try to run that code I get an error message. – Richard Erickson Oct 30 '15 at 14:26
• @RichardErickson. Yep, that was a typo. It should have been n=100 or some such. I just edited it and one other typo/misrendered character, so that the code above now works. – Josh O'Brien Oct 30 '15 at 15:22
• Sorry about the typos. I took it directly from Simon Wood's example but didn't try it out myself because it wasn't exactly what I wanted. – RNB Nov 2 '15 at 5:41

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

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

• Thanks. The main difference is that my main exposure/independent variable is continuous and not binary. But I think I could take a subset of the dataset that just contains two values of the exposure. For example, dat2 <- subset(dat1, age == 20 | age == 25). Then from there create a dichotomous factor variable to compare them. Do you think that will work? – RNB Nov 2 '15 at 5:48
• I just tried it out, but realized that you get a completely different risk ratio than you should from the model, because it relies on crude, unadjusted prevalences. I actually need to figure out how to do this from the model predicted probabilities. Thanks though! – RNB Nov 2 '15 at 7:41