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: