# Comparison of two regression types on the same set of models

I wrote a relatively simple R script to run a bunch of regressions. To be precise, it is a hurdle negative binomial regression (fancy stuff for count data). To print it out easily to Latex readable format, I store the results in a list.

I run five regressions in a row: one baseline regression, two specifications that add one variable at a time to the baseline, and finally two specifications that add one variable at a time to the four regressions before. The first two and the last two additional variables have the same beginning ("A" and "B" resp.) while the variables of the additional specifications have the same ending ("xyz" and "abc", resp.). To accomplish this, I construct the regression formula in a loop.

I also compare these regression results to the results of an ordinary negative binomial regression using Vuong's test (more fancy stuff on count data) right thereafter.

# Load packages
require(MASS) # for NegBin regression
require(pscl) # for Hurdle NegBin regression

Master <- within(Master, {
nofish <- factor(nofish)
camper <- factor(camper)
A_xyz <- persons
A_abc <- persons*2
B_xyz <- livebait
B_abc <- livebait*2
})

# Define variables that are used throughout the script
status <- "without"
var_control <- c("child", "camper")
var_A <- c("A_xyz", "A_abc")
var_B <- c("B_xyz", "B_abc")

estimation_hurdle <- estimation_negbin <- list()

# Define a file to write out the results of Vuong's test
sink(file=paste0("filepath/", status, "/filename.txt"))

# Run baseline regression
writeLines("baseline")
estimation_hurdle$baseline <- hurdle(as.formula(formula_hurdle), data=Master, dist="negbin", link="probit" ) estimation_negbin$baseline <- glm.nb(as.formula(formula_hurdle),
data=Master
)
vuong(estimation_hurdle$baseline, estimation_negbin$baseline)

# Run regression for A variables
for (var in var_A) {
writeLines(paste0("\n\n", var))
estimation_hurdle[[var]] <- hurdle(as.formula(paste(formula_hurdle, var, sep=" + ")),
)
estimation_negbin[[var]] <- glm.nb(as.formula(paste(formula_hurdle, var, sep=" + ")),
data=Master
)
vuong(estimation_hurdle[[var]], estimation_negbin[[var]])
}
rm(var)

# Run regression for B and corresponding A variables
for (i in 1:length(var_B)) {
writeLines(paste0("\n\n", var_B[i]))
estimation_hurdle[[var_B[i]]] <- hurdle(as.formula(paste(formula_hurdle, var_A[i], var_B[i], sep=" + ")),
)
estimation_negbin[[var_B[i]]] <- hurdle(as.formula(paste(formula_hurdle, var_A[i], var_B[i], sep=" + ")),
data=Master
)
vuong(estimation_hurdle[[var_B[i]]], estimation_negbin[[var_B[i]]])
}
rm(i)
sink()

# Additional stuff to write out the results which are not much interesting here


Any idea to improve this code and make it easier/understandable?