# Multinomial tests of heterogeneity using robust standard errors

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 sandwich and 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 factor2:site interactions.

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

1. 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?

2. Should the model based and robust inference be roughly the same for all of these models?

3. Is it okay to say in any of the above cases that what was fit was multinomial logistic regression with robust standard errors?