# \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 <- 0.5
a2 <- 0
a3 <- -1
a4 <- 2
a5 <- 1
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)
}
}
Y <- rpois(n, mu)
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:
pois_mod <- poisRegMisrepEM(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(pois_mod)
# Coefficients:
# Estimate Std. Error z value Pr(>|z|)
# (Intercept) 1.03519 0.02238 46.25615 <2e-16 ***
# X1 0.49875 0.01297 38.45157 <2e-16 ***
# X2b -0.00007 0.01324 -0.00500 0.99601
# X2c -0.98438 0.01926 -51.10084 <2e-16 ***
# X3 1.97794 0.00878 225.20267 <2e-16 ***
# V_star 0.99484 0.01290 77.14885 <2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ---
# AIC AICc BIC
# 4170.836 4170.949 4205.190
# ---
# Log-Likelihood
# -2078.418
# ---
# Lambda: 0.1039615 std.err: 0.01613403
# Fitting an interaction between X2 and X3;
a6 <- -0.5
a7 <- -0.5
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])
}
}
}
Y <- rpois(n, mu)
data$Y <- Y
pois_mod <- poisRegMisrepEM(formula = Y ~ X1 + X2 + X3 + V_star + X2*X3,
v_star = "V_star", data = data)
summary(pois_mod)
# Coefficients:
# Estimate Std. Error z value Pr(>|z|)
# (Intercept) 0.98723 0.02917 33.84255 <2e-16 ***
# X1 0.50135 0.01540 32.56094 <2e-16 ***
# X2b -0.03643 0.03655 -0.99648 0.31902
# X2c -1.02315 0.05170 -19.79103 <2e-16 ***
# X3 1.99527 0.01319 151.22592 <2e-16 ***
# V_star 1.00917 0.01531 65.93335 <2e-16 ***
# X2b:X3 -0.47260 0.02137 -22.11569 <2e-16 ***
# X2c:X3 -0.49639 0.03018 -16.44530 <2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ---
# AIC AICc BIC
# 4096.533 4096.714 4140.702
# ---
# Log-Likelihood
# -2039.266
# ---
# Lambda: 0.1072814 std.err: 0.0162925
# Model fitting with a polynomial effect;
a8 <- -1
mu <- mu*exp(a8*X3^2)
Y <- rpois(n, mu)
data$Y <- Y
pois_mod <- poisRegMisrepEM(formula = Y ~ X1 + X2 + X3 + V_star + X2*X3 + I(X3^2),
v_star = "V_star", data = data)
summary(pois_mod)
# Coefficients:
# Estimate Std. Error z value Pr(>|z|)
# (Intercept) 1.03291 0.04647 22.22701 <2e-16 ***
# X1 0.43783 0.03453 12.68058 <2e-16 ***
# X2b -0.08042 0.05600 -1.43609 0.15098
# X2c -1.02676 0.07523 -13.64912 <2e-16 ***
# X3 2.03183 0.06317 32.16597 <2e-16 ***
# V_star 0.98563 0.03415 28.86175 <2e-16 ***
# I(X3^2) -0.99795 0.03529 -28.27715 <2e-16 ***
# X2b:X3 -0.45828 0.06499 -7.05189 <2e-16 ***
# X2c:X3 -0.47648 0.08912 -5.34623 <2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ---
# AIC AICc BIC
# 3269.698 3269.920 3318.775
# ---
# Log-Likelihood
# -1624.849
# ---
# Lambda: 0.108672 std.err: 0.02181499
# }
Run the code above in your browser using DataLab