# Simulate data
n <- 2000
p0 <- 0.25
X1 <- rbinom(n, 1, 0.4)
X2 <- rnorm(n, 0, 1)
X3 <- rbeta(n, 2, 1)
theta0 <- 0.3
V <- rbinom(n,1,theta0)
V_star <- V
V_star[V==1] <- rbinom(sum(V==1),1,1-p0)
a0 <- 1
a1 <- 2
a2 <- 0
a3 <- 4
a4 <- 2
mu <- exp(a0 + a1*X1 + a2*X2 + a3*X3 + a4*V)
phi <- 0.2
alpha0 <- 1/phi
beta <- 1/mu/phi
Y <- rgamma(n, alpha0, beta)
data <- data.frame(Y = Y, X1 = X1, X2 = X2, X3 = X3, V_star = V_star)
# Split data into training and testing sets
train <- data[1:1800,]
test <- data[1801:2000,]
gamma_fit <- gammaRegMisrepEM(formula = Y ~ X1 + X2 + X3 + V_star,
v_star = "V_star", data = train)
# Predict on test set;
preds <- predict(gamma_fit, newdata = test)
Run the code above in your browser using DataLab