# \donttest{
set.seed(314159)
# Simulate data
n <- 1000
p0 <- 0.25
X1 <- rbinom(n, 1, 0.4)
X2 <- sample(x = c("a", "b", "c"), size = n, replace = TRUE)
X3 <- rnorm(n, 0, 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 <- -1
a4 <- 4
a5 <- 2
mu <- rep(0, n)
for(i in 1:n){
mu[i] <- exp(a0 + a1*X1 + a4*X3 + a5*V )[i]
if(X2[i] == "a" || X2[i] == "b"){
mu[i] <- mu[i]*exp(a2)
}else{
mu[i] <- mu[i]*exp(a3)
}
}
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)
# "a" is the reference
data$X2 <- as.factor(data$X2)
# Model with main effects:
gamma_mod <- gammaRegMisrepEM(formula = Y ~ X1 + X2 + X3 + V_star,
v_star = "V_star", data = data)
# The prevalence of misrepresentation;
(theta0 * p0) / (1 - theta0*(1-p0)) # 0.09677419
# Parameter estimates and estimated prevalence of
# misrepresentation (lambda);
summary(gamma_mod)
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 0.99356 0.03013 32.97245 <2e-16 ***
# X1 2.02152 0.03078 65.68276 <2e-16 ***
# X2b -0.00679 0.03708 -0.18309 0.85477
# X2c -1.02578 0.03684 -27.84599 <2e-16 ***
# X3 3.97883 0.01495 266.21973 <2e-16 ***
# V_star 2.00437 0.03107 64.51234 <2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ---
# AIC AICc BIC
# 5650.696 5650.841 5689.958
# ---
# Log-Likelihood
# -2817.348
# ---
# Lambda: 0.1083894 std.err: 0.01160662
# Fitting an interaction between X2 and X3;
a6 <- -2
a7 <- 2
for(i in 1:n){
if(X2[i] == "c"){
mu[i] <- mu[i]*exp(a6*X3[i])
}else{
if(X2[i] =="b"){
mu[i] <- mu[i]*exp(a7*X3[i])
}
}
}
beta <- 1/mu/phi
Y <- rgamma(n, alpha0, beta)
data$Y <- Y
gamma_mod <- gammaRegMisrepEM(formula = Y ~ X1 + X2 + X3 + V_star + X2*X3,
v_star = "V_star", data = data)
summary(gamma_mod)
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 0.96205 0.03086 31.17145 <2e-16 ***
# X1 2.00411 0.03061 65.46734 <2e-16 ***
# X2b -0.00987 0.03682 -0.26818 0.78862
# X2c -0.99957 0.03733 -26.77449 <2e-16 ***
# X3 3.98282 0.02484 160.31083 <2e-16 ***
# V_star 2.01107 0.03077 65.36550 <2e-16 ***
# X2b:X3 1.95884 0.03573 54.82466 <2e-16 ***
# X2c:X3 -1.98595 0.03567 -55.67827 <2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ---
# AIC AICc BIC
# 5633.984 5634.207 5683.062
# ---
# Log-Likelihood
# -2806.992
# ---
# Lambda: 0.1131951 std.err: 0.01181678
# Model fitting with a polynomial effect;
a8 <- -0.5
mu <- mu*exp(a8*X3^2)
beta <- 1/mu/phi
Y <- rgamma(n, alpha0, beta)
data$Y <- Y
gamma_mod <- gammaRegMisrepEM(formula = Y ~ X1 + X2 + X3 + V_star + X2*X3 + I(X3^2),
v_star = "V_star", data = data)
summary(gamma_mod)
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 1.04312 0.03164 32.96624 <2e-16 ***
# X1 2.04411 0.02929 69.79020 <2e-16 ***
# X2b -0.10418 0.03512 -2.96620 0.00309 **
# X2c -1.08910 0.03531 -30.84683 <2e-16 ***
# X3 4.00265 0.02421 165.31001 <2e-16 ***
# V_star 1.98741 0.02951 67.35719 <2e-16 ***
# I(X3^2) -0.51152 0.01350 -37.90112 <2e-16 ***
# X2b:X3 1.98709 0.03598 55.22750 <2e-16 ***
# X2c:X3 -2.03395 0.03692 -55.09491 <2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ---
# AIC AICc BIC
# 4559.969 4560.236 4613.954
# ---
# Log-Likelihood
# -2268.984
# ---
# Lambda: 0.111464 std.err: 0.01173143
# }
Run the code above in your browser using DataLab