I am looking at differences in frequencies of categorical variables collected at different sites, adjusting for stratification factors of age and sex. I would like to use robust standard error estimation as these communities are small and there are unmeasured sources of correlation among participants, but the sample is population representative. The sample size is relatively large. Age is measured continuously. A simulated dataset could be:
set.seed(123) n <- 300 age <- rnorm(n, 60, 10) sex <- sample(c('male', 'female'), n, replace=T) site <- sample(0:1, n, replace=T) factor1 <- sample(0:1, n, replace=T) factor2 <- sample(0:2, n, replace=T)
For binary factors, the obvious logistic regression model is easily calculated and robust error estimation is implemented for these models. This is easily done in R using the
lmtest packages. For instance:
library(sandwich) library(lmtest) vhc0 <- function(x) vcovHC(x, type='HC0') ## no correction logreg <- glm(factor1 ~ site + age + sex, family=binomial) coeftest(logreg, vcov=vhc0)['site', ] > coeftest(logreg, vcov=vhc0)['site', ] Estimate Std. Error z value Pr(>|z|) -0.1831499 0.2321805 -0.7888254 0.4302141
A nearly equivalent robust and model based inference is obtained using log-linear models via:
Freq <- as.data.frame(table(factor1, site, age, sex))) Freq$age <- as.numeric(Freq$age) loglinreg <- glm(Freq ~ factor1*(site + age + sex), data=Freq, family=poisson) coeftest(loglinreg, vcov=vhc0)['factor11:site1', ] > coeftest(loglinreg, vcov=vhc0)['factor11:site1', ] Estimate Std. Error z value Pr(>|z|) -0.1868107 0.2162724 -0.8637748 0.3877116
Which I theoretically believe should be asymptotically equivalent.
Now the question of inference about
factor2 is not so clear cut. Without adjustment and without robust standard errors, the chisquare test of independence is the obvious choice. With adjustment, a multinomial model could be used for an equivalent. But no multinomial models in R provide robust error estimation, except using by-hand implementations of the bootstrap.
Option 1) fit the plain vanilla logistic regression, except by doing sequential analysis comparing each factor level of
factor2 to the referent value e.g. fitting the 2:1 model and 3:1 model, and taking the minimum of their respective $p$-values for inference about
site to evaluate the overall test of heterogeneity... but this will still be of relatively low power.
fit1.0 <- glm(factor2 == 1 ~ site + age + sex, family=binomial, subset=factor2 != 2) fit2.0 <- glm(factor2 == 2 ~ site + age + sex, family=binomial, subset=factor2 != 1) c(coef(fit1.0)['site'], coef(fit2.0)['site']) min( coeftest(fit1.0, vcov=vhc0)['site', 'Pr(>|z|)'], coeftest(fit2.0, vcov=vhc0)['site', 'Pr(>|z|)'] ) site site 0.02480113 0.40289471  0.1694746
Option 2) fit the log linear model in a manner similar to the above and performing joint inference on the
Freq <- as.data.frame(table(factor2, site, age, sex)) Freq$age <- as.numeric(Freq$age) loglinreg <- glm(Freq ~ factor2*(site + age + sex), data=Freq, family=poisson) loglinreg.0 <- glm(Freq ~ factor2*(age + sex), data=Freq, family=poisson) coef(loglinreg)[paste0('factor2', 1:2, ':site1')] waldtest(loglinreg, loglinreg.0, vcov=vhc0) factor21:site1 factor22:site1 -0.003194891 0.391984947 Wald test Model 1: Freq ~ factor2 * (site + age + sex) Model 2: Freq ~ factor2 * (age + sex) Res.Df Df F Pr(>F) 1 3588 2 3591 -3 0.922 0.4292
Option 3) fit a "stacked" logistic regression model which stratifies by multinomial comparison
site.2 <- rep(site, 2) age.2 <- rep(age, 2) sex.2 <- rep(sex, 2) factor2.10 <- factor2 == 1 factor2.20 <- factor2 == 2 factor2.10[factor2 == 2] <- NA factor2.20[factor2 == 1] <- NA factor2.2 <- c(factor2.10, factor2.20) analysis <- factor(rep(c('1:0', '2:0'), each=n)) lrstack <- glm(factor2.2 ~ analysis*(age.2 + sex.2 + site.2), family=binomial) lrstack.0 <- glm(factor2.2 ~ analysis*(age.2 + sex.2 ), family=binomial) coef(lrstack)[c('site.2', 'analysis2:0:site.2')] waldtest(lrstack, lrstack.0, vcov=vhc0) site.2 analysis2:0:site.2 0.02480113 0.37809359 Wald test Model 1: factor2.2 ~ analysis * (age.2 + sex.2 + site.2) Model 2: factor2.2 ~ analysis * (age.2 + sex.2) Res.Df Df F Pr(>F) 1 386 2 388 -2 0.9475 0.3886
Should the model based coefficients from a log linear model, a stacked logistic regression, sequential logistic regression, and a multinomial model be the same when adjusting for other covariates?
Should the model based and robust inference be roughly the same for all of these models?
Is it okay to say in any of the above cases that what was fit was multinomial logistic regression with robust standard errors?