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:


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

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

  • $\begingroup$ 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. $\endgroup$ – Richard Erickson Oct 30 '15 at 14:26
  • $\begingroup$ @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. $\endgroup$ – Josh O'Brien Oct 30 '15 at 15:22
  • $\begingroup$ 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. $\endgroup$ – 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

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"))
## 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")
| cite | improve this answer | |
  • $\begingroup$ 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? $\endgroup$ – RNB Nov 2 '15 at 5:48
  • $\begingroup$ 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! $\endgroup$ – RNB Nov 2 '15 at 7:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.