# 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)
gamma_fit <- gammaRegMisrepEM(formula = Y ~ X1 + X2 + X3 + V_star,
v_star = "V_star", data = data)
summary(gamma_fit)
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 1.00137 0.03413 29.33857 <2e-16 ***
# X1 2.01388 0.02154 93.48440 <2e-16 ***
# X2 -0.00193 0.01038 -0.18589 0.85255
# X3 4.00101 0.04560 87.74528 <2e-16 ***
# V_star 2.00567 0.02240 89.54515 <2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ---
# AIC AICc BIC
# 23362.50 23362.56 23401.71
# ---
# Log-Likelihood
# -11674.25
# ---
# Lambda: 0.09635239 std.err: 0.007641834
Run the code above in your browser using DataLab