I need to report the mean, and SD of two groups for a specific variable, as well as the odds ratio of the two groups using robust logistic regression. The variable does have outliers. I instead want to report the robust mean and sd with a method that is similarly used in the robust logistic regression.
I simply want to get the new robust calculated mean and SD of the first term using the same methods from robustbase::glmrob R function. Method is default (method="Mqle") and family=binomial(link = "logit"). I calculated the Huber M-estimator of the mean and SD using MASS::hubers(x, k=1.3450) but I am not sure this is the estimator method used in robustbase::glmrob. Do you know what method is used and how I can use that method to calculate the robust mean and SD in R? Based on reading the original papers, it seems that to get a robust estimator of variance, they choose a simple M-estimator (Huber’s proposal 2). This variance estimator is calculated using MASS::hubers(), not MASS::huber(). The latter is based on the median absolute deviation (MAD) scale. However, glmrob does not seem to use a robust mean in the formulas. I still believe its better to report robust means and SD so that the summary stats and models are more comparable.
method="Mqle" fits a generalized linear model using Mallows (Mallows quasi-likelihood estimator) or Huber type robust estimators, as described in Cantoni and Ronchetti (2001) and Cantoni and Ronchetti (2006).
Eva Cantoni and Elvezio Ronchetti (2001) Robust Inference for Generalized Linear Models. JASA 96 (455), 1022–1030.
Eva Cantoni and Elvezio Ronchetti (2006) A robust approach for skewed and heavy-tailed outcomes in the analysis of health care expenditures. Journal of Health Economics 25, 198–213.
Here is an example:
Here is an example dataframe with a binomial response variable "Group", a numeric variable and covariable.
> df <- data.frame(group = as.factor(c(rep(1, 10), rep(2,10))),
variable = c(2,6,1,2,2,2,4,0,2,6,7,19,5,6,0,10,5,2,5,0),
covariable = c(19,22,16,18,18,19,19,15,18,13,16,16,19,17, 21,12,14,19,18,18)
)
> df
group variable covariable
1 1 2 19
2 1 6 22
3 1 1 16
4 1 2 18
5 1 2 18
6 1 2 19
7 1 4 19
8 1 0 15
9 1 2 18
10 1 6 13
11 2 7 16
12 2 19 16
13 2 5 19
14 2 6 17
15 2 0 21
16 2 10 12
17 2 5 14
18 2 2 19
19 2 5 18
20 2 0 18
Calculate the standard mean and SD
> aggregate(variable ~ group, df, function(x) c(mean = mean(x), sd = sd(x)))
group variable.mean variable.sd
1 1 2.7 2.0
2 2 5.9 5.6
Using robust logistic regression. Notice that rows 15 and 20 are outliers and are thus weighted. These outliers are due to the variable, not covariable. I know this because the weights hardly changed when I removed the covariable term.
> robustbase::glmrob(group ~ variable + covariable , data=df, family = binomial(link = "logit")) %>% summary()
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.264760 4.218720 -0.300 0.764
variable 0.280474 0.208623 1.344 0.179
covariable 0.006005 0.220648 0.027 0.978
Robustness weights w.r * w.x:
18 weights are ~= 1. The remaining 2 ones are
15 20
0.7611 0.7543
Number of observations: 20
Fitted by method ‘Mqle’ (in 7 iterations)
(Dispersion parameter for binomial family taken to be 1)
No deviance values available
Algorithmic parameters:
acc tcc
0.0001 1.3450
maxit
50
test.acc
"coef"
Calculate the Huber M-estimator of the mean and "SD". Notice that the mean and SDs are different from the original. I "think" that the same M-estimator method used in MASS::hubers(y, k=1.3450) to estimate the mean and variance is used in glmrob(method="Mqle") function. If not, please correct me.
> aggregate(variable ~ group, df, function(y) c(mean.h = MASS::hubers(y, k=1.3450)$mu, sd.h = MASS::hubers(y, k=1.3450)$s))
group variable.mean.h variable.sd.h
1 2.6 2.2
2 5.1 4.4
