QRRT: Poisson regression for QRRT data
In meca7653/QRRT: Poisson regression with QRRT data

Description Usage Arguments Value Author(s) References Examples

View source: R/QRRT_V7.R

Poisson regression for quantitative randomized response technique (QRRT) data, based on maximum likelihood implemented via the EM algorithm

1
2
3

QRRT(Formula, Data, Disperse = 1, beta = NULL, n_times = 1,
  offset = NULL, b_distribution = c(6, 7, 4, 2, 2, 1, 1, 1, 1, 25)/50,
  tolerance = 1e-08)

`Formula`	symbolic description of the model to be fitted using the same specification as the lm function for linear models
`Data`	data frame containing the variables in the model
`Disperse`	standard deviation of the random starting values for beta
`beta`	NULL for random starting values or a specified starting value for regression coefficients
`n_times`	number of random starting values for beta
`offset`	NULL or a numeric vector of length equal to the number of cases
`b_distribution`	probability distribution on 0, 1, ..., m, m + 1. Default is from Conteh et al. (2015).
`tolerance`	EM convergence criterion for change in log-likelihood

List containing the following elements:

`Results`	a matrix including maximum likelihood estimates of regression coefficients and estimated asymptotic standard errors, along with corresponding t-statistics and p-values
`AIC`	Akaike's An Information Criterion: minus twice the maximized log-likelihood plus twice the number of parameters
`Maximized_Log_Likelihood`	log-likelihood at the estimated coefficient values
`Likelihood_by_Starting_Value`	maximized likelihood at each of the n_times starting values
`Number_Missing_Values`	number of missing value in Data
`Score`	score vector evaluated at the maximum likelihood estimate of regression coefficient vector
`Inverse_Fisher_Information`	estimated variance-covariance matrix for maximum likelihood estimate of regression coefficient vector

Meng Cao, F. Jay Breidt

Conteh A, Gavin MC, Solomon J. Quantifying illegal hunting: A novel application of the quantitative randomised response technique. Biological Conservation. 2015;189:16-23.

Liu P, Chow L. A New Discrete Quantitative Randomized Response Model. Journal of the American Statistical Association. 1976;71(353):72<e2><80><93>73.

Cao M, Breidt F.J, Solomon J, Conteh A, Gavin M. Understanding the Drivers of Sensitive Behavior Using Poisson Regression from Quantitative Randomized Response Technique Data. Paper under written. Colorado State University

#-------------------------------------------------------------------------------------
# Simulate a QRRT data set:
#-------------------------------------------------------------------------------------
# True count is Poisson with mean lambda =
# exp(beta_0+beta_1*x1+beta_2*x2+beta_3*dummyB+beta_4*dummyC+beta_5*x1*x2)
# where we have continuous predictors x1 and x2 and a categorical predictor
# x3 with levels A, B, C.  The categorical variable is encoded
# as a dummy variable for level B and a dummy variable for level C,
# with level A as baseline.
beta <- cbind(c(1.5, 1.0, -0.5, 0.4, 0.3, 0.2))
n <- 400 # sample size
set.seed(12345)
x1 <- sample(15:70, size = n, replace = TRUE) / 50
x2 <- rnorm(n, 0, 1)
x3 <- sample(c("A", "B", "C"), size = n, replace = TRUE)
dummyB <- 1 * (x3 == "B")
dummyC <- 1 * (x3 == "C")
X <- cbind(rep(1, n), x1, x2, dummyB, dummyC, x1 * x2)
lambda <- exp( X %*% beta)
Ti <- rpois(n, lambda) # true response T_i
range(Ti)
# Using same known distribution as Conteh et al. (2015).
# This is a distribution on 0, 1, ..., m, m + 1.
b_distribution <- c(6, 7, 4, 2, 2, 1, 1, 1, 1, 25) / 50
sum(b_distribution)
m <- length(b_distribution) - 2
# Draw the random integers B_i ~ b(r) for i = 1, 2, ..., n.
Bi <- sample(0:(m + 1), size = n, replace = TRUE, prob = b_distribution)
# If B_i < m + 1, observe B_i; otherwise, observe truth.
Ri <- Bi * (Bi < (m + 1)) + Ti * (Bi == (m + 1))
Sim_Data <- data.frame(x1, x2, x3, Ri)
# Simulated data are now complete.

#------------------------------------------------------------------------------------
# Now use the QRRT code to fit models to the simulated data.
#------------------------------------------------------------------------------------
# Fit model with all two-way interactions, using 10 random starting values.
fit_2way <- QRRT(Formula = Ri ~ (x1 + x2 + as.factor(x3)) ^ 2,
                 Data = Sim_Data, Disperse = 1, beta = NULL, n_times = 10,
                 offset = NULL,
                 b_distribution = c(6, 7, 4, 2, 2, 1, 1, 1, 1, 25) / 50)
fit_2way$Results
#------------------------------------------------------------------------------------
#Fit model correctly-specified as the actual simulation model.
fit_truemodel <- QRRT(Formula = Ri ~ (x1 + x2) ^ 2 + as.factor(x3),
                      Data = Sim_Data, Disperse = 1, beta = NULL, n_times = 10,
                   offset = NULL,
                   b_distribution = c(6, 7, 4, 2, 2, 1, 1, 1, 1, 25) / 50)
fit_truemodel$Results
#------------------------------------------------------------------------------------
# Coefficient of x1 in simulated model is 1; use this fact to illustrate
# inclusion of an offset in the model.
fit_offset <- QRRT(Formula = Ri ~ x2 + x1:x2 + as.factor(x3),
                Data = Sim_Data, Disperse = 1, beta = NULL, n_times = 10,
                offset = x1,
                b_distribution = c(6, 7, 4, 2, 2, 1, 1, 1, 1, 25) / 50)
fit_offset$Results
#------------------------------------------------------------------------------------
round(fit_2way$Results, 3)

#------------------------------------------------------------------------------------
# Real data example
data(QRRT_data)
F.formula <- times.month.trap ~ Residencetime + Education.level + HouseholdSize +
 PAAfterWar + NoRestriction + EfficientConservation + AfterPatrol +
 QuickApprehend + Punished + OutsidersHunted + ResidentsHunted +
 AnimalsRare + Seaside + Rural + HighDisplace
res <- QRRT(Formula = F.formula,
 Data = QRRT_data,
 n_times = 1)