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?
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) :)