Rank-deficient models

[rvlenth]

Classification:

Enhancement suggestion

Summary

An emmeans user posted a question about emmprep whereby he wasn't able to use emmprep() for a model because it had nested fixed effects, and that caused rank deficiencies in the model.

Reproducible Example (if applicable)

library(metafor)
## Loading required package: Matrix
## Loading required package: metadat
## Loading required package: numDeriv
## 
## Loading the 'metafor' package (version 4.2-0). For an
## introduction to the package please type: help(metafor)
library(emmeans)

dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg)

dat <- transform(dat,
                 acat = factor(ablat < 40, labels = c("high", "low")),
                 era = factor(1+as.integer((year-1940)/15), 
                              labels = c("early","mid","late")))

with(dat, table(acat, era))
##       era
## acat   early mid late
##   high     3   2    1
##   low      0   2    5

mod <- rma(yi, vi, mods = ~ acat * era, data = dat)
## Warning: Redundant predictors dropped from the model.
mod
## 
## Mixed-Effects Model (k = 13; tau^2 estimator: REML)
## 
## tau^2 (estimated amount of residual heterogeneity):     0.0854 (SE = 0.0810)
## tau (square root of estimated tau^2 value):             0.2923
## I^2 (residual heterogeneity / unaccounted variability): 68.48%
## H^2 (unaccounted variability / sampling variability):   3.17
## R^2 (amount of heterogeneity accounted for):            72.73%
## 
## Test for Residual Heterogeneity:
## QE(df = 8) = 23.1259, p-val = 0.0032
## 
## Test of Moderators (coefficients 2:5):
## QM(df = 4) = 18.5649, p-val = 0.0010
## 
## Model Results:
## 
##                 estimate      se     zval    pval    ci.lb    ci.ub      
## intrcpt          -0.9706  0.2437  -3.9823  <.0001  -1.4483  -0.4929  *** 
## acatlow           1.1731  0.3627   3.2348  0.0012   0.4623   1.8839   ** 
## eramid           -0.3952  0.4240  -0.9320  0.3513  -1.2262   0.4359      
## eralate          -0.4710  0.4060  -1.1600  0.2461  -1.2667   0.3248      
## acatlow:eramid   -0.1059  0.6060  -0.1747  0.8613  -1.2937   1.0819      
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

RG <- emmprep(mod)
## Error in emmprep(mod): Cannot use function when some redundant predictors were dropped from the model.

# compare with lm:
LM <- lm(yi ~ acat * era, data = dat)
coef(LM)
##     (Intercept)         acatlow          eramid         eralate  acatlow:eramid 
##      -1.0869385       1.0049180      -0.2727705      -0.3546127       0.3430389 
## acatlow:eralate 
##              NA

(RG <- ref_grid(LM))
## 'emmGrid' object with variables:
##     acat = high, low
##     era = early, mid, late
## Some things are non-estimable (null space dim = 1)

confint(RG)
##  acat era   prediction    SE df lower.CL upper.CL
##  high early    -1.0869 0.331  8   -1.851   -0.323
##  low  early     nonEst    NA NA       NA       NA
##  high mid      -1.3597 0.406  8   -2.295   -0.424
##  low  mid      -0.0118 0.406  8   -0.947    0.924
##  high late     -1.4416 0.574  8   -2.765   -0.118
##  low  late     -0.4366 0.257  8   -1.028    0.155
## 
## Confidence level used: 0.95

^{Created on 2023-07-22 with reprex v2.0.2}

^{Created on 2023-07-22 with reprex v2.0.2}

Notes

The emmeans package is capable of supporting rank-deficient models. However, it needs information about which predictors were dropped. That information is also useful to some users. With some small changes in the object created it would be possible to enable that support and take advantage of emmeans's ability to handle those cases rather than just issuing an error message.

The conventional way to handle rank deficiencies is the way lm() does it. The cofficients for all the predictors are returned by coef(), and the ones that are dropped are given the value NA. The vcov() function, on the other hand, just returns the covariances for the included predictors. Note in the confidence intervals for RG (the reference grid for LM), we were able to identify the non-estimable case (the cell that was empty). It would seem that if rma() also returned the NAs as coefficients for the dropped predictors, that is all it would take, as the qdrg() function can figure this stuff out too).

One point of common confusion (or at least to me, a number of years ago): The NA coefficients do not actually represent missing values! Instead, they represent coefficients that were set to zero. This is clear as soon as you make the connection that dropping a predictor is the same thing as setting its coefficient to zero. When there are rank deficiencies, the solution to the normal equations is not unique; and in R, the particular solution returned is the one with those constraints to zero. There are many other valid solutions, and most of them give non-zero weight to those predictors. So we really are not saying that the dropped predictors are unimportant -- just that we didn't use them for this one particular solution.

sessionInfo()

R version 4.3.1 (2023-06-16 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 11 x64 (build 22621)

Matrix products: default


locale:
[1] LC_COLLATE=English_United States.utf8  LC_CTYPE=English_United States.utf8   
[3] LC_MONETARY=English_United States.utf8 LC_NUMERIC=C                          
[5] LC_TIME=English_United States.utf8    

time zone: America/Chicago
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] metafor_4.2-0        numDeriv_2016.8-1.1  metadat_1.2-0        Matrix_1.5-4.1      
[5] emmeans_1.8.7-090005

loaded via a namespace (and not attached):
 [1] styler_1.10.0      sandwich_3.0-2     utf8_1.2.3         generics_0.1.3     lattice_0.21-8    
 [6] digest_0.6.33      magrittr_2.0.3     evaluate_0.21      grid_4.3.1         estimability_1.4.1
[11] mvtnorm_1.1-3      fastmap_1.1.1      R.oo_1.25.0        R.cache_0.16.0     processx_3.8.1    
[16] R.utils_2.12.2     survival_3.5-5     ps_1.7.5           multcomp_1.4-23    purrr_1.0.1       
[21] fansi_1.0.4        scales_1.2.1       TH.data_1.1-2      codetools_0.2-19   cli_3.6.1         
[26] rlang_1.1.1        R.methodsS3_1.8.2  reprex_2.0.2       munsell_0.5.0      splines_4.3.1     
[31] yaml_2.3.7         withr_2.5.0        tools_4.3.1        parallel_4.3.1     coda_0.19-4       
[36] dplyr_1.1.2        colorspace_2.1-0   ggplot2_3.4.2      mathjaxr_1.6-0     vctrs_0.6.3       
[41] R6_2.5.1           zoo_1.8-12         lifecycle_1.0.3    fs_1.6.2           MASS_7.3-60       
[46] callr_3.7.3        clipr_0.8.0        pkgconfig_2.0.3    pillar_1.9.0       gtable_0.3.3      
[51] glue_1.6.2         xfun_0.39          tibble_3.2.1       tidyselect_1.2.0   highr_0.10        
[56] rstudioapi_0.14    knitr_1.43         farver_2.1.1       xtable_1.8-4       htmltools_0.5.5   
[61] nlme_3.1-162       rmarkdown_2.21     labeling_0.4.2     compiler_4.3.1  

[rvlenth]

Follow-up: It turns out I already had code to figure out where to fill-in the NAs. I think you can just take the check on coef.na and it'll work just fine. See the answer I posted on the SO site I referenced above.

PS -- there was a glitch in the code for qdrg that made it not always get the V matrix right in rank-deficient models, but I don't think it affects your models as it only comes into play when the number of rows/columns is too large. But I will be sending an update within a week or two that will have that fixed.

[wviechtb]

Hi Russell -- thanks for the suggestion. I am just coming back to the office from a 5-week absence due to conferences and holidays. It will take me a while to get back to this, but will respond when I find the time.