amrei-stammann / alpaca

An R-package for fitting glm's with high-dimensional k-way fixed effects
43 stars 6 forks source link

getAPEs and clustered SEs #20

Closed eohne closed 2 months ago

eohne commented 3 months ago

Hi, first thanks for the awesome package. I have a quick question. Does the summary output of getAPEs actually allow for clustering?

I run a Logit model with two fixed effects using feglm.

res = feglm(Y ~ x1 + x1 |fe1 + fe2| cluster_var + other_cluster_var,family=binomial(link = "logit"))

I next use the biasCorr command.

res = biasCorr(res)

If I use summary(bias_corr_res, cluster = ~cluster_var, type="clustered") I get clustered standard errors back and all is okay. i.e.

summary(bias_corr_res, cluster = ~cluster_var, type="clustered")`

If I instead want the output for getAPEs instead standard errors are not clustered. The command I run is this:

summary(getAPEs(bias_corr_res), cluster=~cluster_var, type="clustered")

which gives the same standard error as:

summary(getAPEs(bias_corr_res), cluster=~other_cluster_var, type="clustered")

or this:

summary(getAPEs(bias_corr_res))

Edit: Another quick question. Is there a way to render the summary table as a LaTeX table? Like etable, stargazer do?

grantmcdermott commented 2 months ago

Another quick question. Is there a way to render the summary table as a LaTeX table? Like etable, stargazer do?

You can write some quick custom tidiers, which would enable integration with the (fantastic) modelsummary package.

library(alpaca)
data(psid, package = "bife")

mod = feglm(
  LFP ~ KID1 + KID2 + KID3 + log(INCH) | ID + TIME,
  data   = psid,
  family = binomial("probit")
  )

apes = getAPEs(mod)
bcapes = getAPEs(biasCorr(mod))

# add custom functions to extract estimates (tidy) and goodness-of-fit (glance) information
tidy.APEs = function(x, ...) {
  cm = summary(x, ...)$cm
  cm = cbind(rownames(cm), data.frame(cm, row.names=NULL))
  ci = data.frame(confint(x, ...), row.names=NULL)
  ret = setNames(
    cbind(cm, ci),
    c('term', 'estimate', 'std.error', 'statistic', 'p.value', 'conf.low', 'conf.high')
  )
  ret
}

glance.APEs = function(x, ...) {
    data.frame(class = 'APE')
}

library(modelsummary)

 modelsummary(
  list('Naive' = apes, 'Bias-corrected' = bcapes),
  output = 'markdown' # change to e.g. 'regtab.tex' for latex writing
)
Naive Bias-corrected
KID1 -0.088 -0.097
(0.008) (0.008)
KID2 -0.045 -0.049
(0.007) (0.007)
KID3 -0.001 -0.001
(0.005) (0.005)
log(INCH) -0.030 -0.034
(0.008) (0.008)
class APE APE

Created on 2024-09-20 with reprex v2.1.1

amrei-stammann commented 2 months ago

Hi,

thanks for your question.

Currently, alpaca does not support clustered SEs for APEs but only for the coefficients. This extension is on our agenda, but requires some careful thoughts and derivations first due to some complications with bias-correted APEs.

I admit that it’s a bit unfortunate that there is no warning, when you are trying to compute clustered SEs for APEs. So it’s good that we aware of it now.

Although not perfect, in the meanwhile it might be ok to test the signicance of the effect of your variable of interest with the bias corrected coefficients and clustered SEs.

Best wishes, Amrei

eohne commented 2 months ago
  1. Thanks @grantmcdermott for suggesting the modelsummary package and the code snippet.
  2. @amrei-stammann I agree, a warning or simply a sentence mentioning this in the docs would be good in the meantime.

Thanks again for your awsome package!

Will close this issue as I now know all that I need to know. Thanks again to you both!

Best wishes Elias