Search code examples
rglm

r betareg predict(...) equivalent to the GLM predict(....)


Context: Using the betareg package to reproduce this example

I think it boils down to getting betareg equivalent of the GLM predict(...)

# Generate data
mydata <- data.frame(Ft = c(1, 6, 11, 16, 21, 2, 7, 12, 17, 22, 3, 8, 
                            13, 18, 23, 4, 9, 14, 19, 5, 10, 15, 20),
                     Temp = c(66, 72, 70, 75, 75, 70, 73, 78, 70, 76, 69, 70, 
                              67, 81, 58, 68, 57, 53, 76, 67, 63, 67, 79),
                     TD = c(0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 
                            0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 0))

# Run logistic regression model
model <- glm(TD ~ Temp, data=mydata, family=binomial(link="logit"))

# Create a temporary data frame of hypothetical values
temp.data <- data.frame(Temp = seq(53, 81, 0.5))

# Predict the fitted values given the model and hypothetical data
predicted.data <- as.data.frame(predict(model, newdata = temp.data, 
                                        type="link", se=TRUE))

glimpse(predicted.data)

  Observations: 57
  Variables: 3
  $ fit            <dbl> 2.73827620, 2.62219483, 2.50611346, 2.39003209, 2.27395072, 2.15786934, 2.0...
  $ se.fit         <dbl> 1.7132157, 1.6620929, 1.6111659, 1.5604536, 1.5099778, 1.4597631, 1.4098372...
  $ residual.scale <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ...

AFAICT the betareg function predict(...) does not produce the se.fit values which are required to display the Confidence Interval ribbon in the chart.

Does anyone know how to mimic the GLM predict(...) in this example?


Solution

  • I have used the package 'effects' and the function allEffects()

    BetaReg <- betareg(Value~x, data = Data)
    Effects <- as.data.frame(allEffects(BetaReg, xlevels=list(x=BetaReg$x)))$x
    Effects
    

    I think it is correct (looks about right on a graph) :)