2

I am trying to bootstrap fitted values of my betareg model in R. I have read many other questions and websites such as https://stats.stackexchange.com/questions/234254/confidence-intervals-for-beta-regression/234256#234256, How to Bootstrap Predictions and Levels of Confidence for Beta Regression Model in R, https://stats.stackexchange.com/questions/86432/how-do-i-predict-with-standard-errors-using-betareg-package-in-r. However, none of these have provided me with an answer I could use for my data. Apparently, calculating the CI for the fitted values of a betareg model is not as straight forward as I thought it would be. I understand that I can use the boot and boot.ci functions from the boot package, but I don't really understand how I should write the statistics function and how to incorporate it in the boot.ci function. Additionally, I have tried the confint function from the betaboost package, but that only gives the 95% CI of the mean, where I am trying to find the CI for my fitted values, so that I can plot the CI together with the model. I'm hoping somebody could show me how to use the bootstrap method to find the 95% CI of the fitted values. Help is much appreciated!

I am investigating the influence of X on Y, both proportions. The data + model looks like this.

graph

My R script

library(dplyr)
library(ggplot2)
library(betareg)
rm(list = ls())

df <- data.frame(propX = c(0.7, 0.671, 0.6795, 0.79, 0.62, 0.62, 0.6413, 0.089, 0.4603, 0.04, 0.0418, 0.46, 0.5995, 0.532, 0.65, 0.6545, 0.74, 0.74, 0.02, 0.02, 0, 0, 0, 0.45, 0.8975, 0.92, 0.898, 0.89, 0.86, 0.69, 0.755, 0.775, 0.585, 0.585, 0.55),
                 propY = c(0.666666666666667, 0.40343347639485, 0.7, 0, 0, 0.0454545454545455, 0.25, 0.707070707070707, 0.629213483146067, 0.882352941176471, 0.942857142857143, 0.451612903225806, 0.0350877192982456, 0.5, 0.484375, 0, 0.0208333333333333, 0.240740740740741, 0.804568527918782, 0.666666666666667, 1, 1, 1, 0.552238805970149, 0.2, 0, 0, 0, 0, 0, 0.12972972972973, 0.0894117647058824, 0.576158940397351, 0, 0),
                 pointWeight = c(3,233,10,89,4,22,44,99,89,17,35,341,57,36,128,39,144,54,394,12,46,229,55,67,5,28,2,160,124,294,555,425,302,116,48))

df$propY <- (((df$propY*(length(df$propY)-1))+0.5)/length(df$propY)) # Transform the data so all data is (0,1)
mybetareg <- betareg(propY ~ propX, data = df, weights = pointWeight, link = "logit")
minoc <- min(df$propX)
maxoc <- max(df$propX)
new.x <- expand.grid(propX = seq(minoc, maxoc, length.out = 1000))
new.y <- predict(mybetareg, newdata = new.x)

# I would like to calculate 95% CI for new.y using the bootstrap method

new.y <- data.frame(new.y)
addThese <- data.frame(new.x, new.y)
addThese <- rename(addThese, propY = new.y)
ggplot(df, aes(x = propX, y = propY)) +
  geom_point(aes(size = pointWeight)) +
  geom_smooth(data = addThese, stat = 'identity') + # here I could then add aes(ymin = lwr, ymax = upr)
  scale_x_continuous(breaks = seq(0,1,0.2), limits = c(0,1)) +
  scale_y_continuous(breaks = seq(0,1,0.2), limits = c(0,1)) +
  theme_bw()
StupidWolf
  • 45,075
  • 17
  • 40
  • 72
Peter
  • 343
  • 5
  • 17

1 Answers1

2

After some trial and error I started working with a different approach to analyse proportional data, namely gam (gam package) with the betar family (mgcv package). This yields exactly the same results as betareg, but it offers more options, such as random effects and standard errors. After the analysis I predict the fitted values and their SE, from which I calculate the 95% confidence interval. The following script should produce a graph with a confidence interval, just fill in your variables and dataset. 

mygam = gam(y ~ x, family=betar(link="logit"), data = df, weights = pointWeight)
min <- min(df$x)
max <- max(df$x)
new.x <- expand.grid(x = seq(min, max, length.out = 1000))
new.y <- predict(mygam, newdata = new.x, se.fit = TRUE, type="response")
new.y <- data.frame(new.y)
addThese <- data.frame(new.x, new.y)
addThese <- rename(addThese, y = fit, SE = se.fit)
addThese <- mutate(addThese, lwr = y - 1.96 * SE, upr = y + 1.96 * SE) # calculating the 95% confidence interval
ggplot(df, aes(x = x, y = y)) +
  geom_point(aes(size = pointWeight)) +
  geom_smooth(data = addThese, aes(ymin = lwr, ymax = upr), stat = 'identity')
Peter
  • 343
  • 5
  • 17