aprstable with robust standard errors for model with singularities
I am trying to create a table of estimation results with robust standard error using the aprstable package. However, because my model has singularities that aren't easily avoided (they come from overlapping factors with hundreds of levels), aprstable refuses to provide the table I am trying to produce. Here's a working example (adapted from a related thread):
library(sandwich) set.seed(101) dat<-data.frame(one=c(sample(1000:1239)), two=c(sample(200:439)), three=c(sample(600:839)), Jan=c(rep(1,20),rep(0,220)), Feb=c(rep(0,20),rep(1,20),rep(0,200)), Mar=c(rep(0,40),rep(1,20),rep(0,180)), Apr=c(rep(0,60),rep(1,20),rep(0,160)), May=c(rep(0,80),rep(1,20),rep(0,140)), Jun=c(rep(0,100),rep(1,20),rep(0,120)), Jul=c(rep(0,120),rep(1,20),rep(0,100)), Aug=c(rep(0,140),rep(1,20),rep(0,80)), Sep=c(rep(0,160),rep(1,20),rep(0,60)), Oct=c(rep(0,180),rep(1,20),rep(0,40)), Nov=c(rep(0,200),rep(1,20),rep(0,20)), Dec=c(rep(0,220),rep(1,20))) model <- lm(one ~ two + three + Jan + Feb + Mar + Apr + May + Jun + Jul + Aug + Sep + Oct + Nov + Dec, data=dat) summary(model) model$se <- vcovHC(model) library('apsrtable') apsrtable(model, se='robust')
executing this code produces
> apsrtable(model, se='robust') Error in s$coefficients[, 3] <- tval <- est/x$se : number of items to replace is not a multiple of replacement length In addition: Warning message: In est/x$se : longer object length is not a multiple of shorter object length
As far as I can tell this error occurs because length(sqrt(diag(model$se))) is 14 while length(coefficients(model)) is 15, i.e. there are more coefficients than there are standard errors on the diagonal of their variance-covariance matrix (see also the relevant bit of source code in the aprstable package).
Is there any easy way to fix this? Insert NA rows and columns into the variance-covariance matrix in the appropriate places before passing it to apsrtable() perhaps? Something more elegant?
This was an easy fix — coef(model) returns the thing with NA, but coef(summary(model)) returns the coef of the summary. I'll submit an update to cran soon, but the package at https://github.com/malecki/apsrtable is updated.