Kleibergen-Paap test code broken; Effective F-statistic implementation

[apoorvalal]

K-P error

Computing the Kleibergen-Paap F-statistic using fitstat(mod, "kpr") produces the following error with the dev version of fixest. The equivalent computation in lfe works fine (and is identical to the documented ivreg2 version in stata).

Error in str2lang(value_raw): <text>:2:0: unexpected end of input
1: (+Intercept
   ^
Traceback:

1. fitstat(m_feols, "kpr")
2. kp_stat(x)
3. proj_on_U(x, Z)
4. resid(feols.fit(Z, U))
5. feols.fit(Z, U)
6. feols(env = env)
7. .xpd(rhs = fml[[3]])
8. value2stringCall(rhs, call = TRUE, check = check)
9. str2lang(value_raw)

MWE:

library(data.table); library(fixest); library(lfe)
set.seed(42)
N = 1e3
df = data.table(
  e = rnorm(N),
  z1 = rbinom(N, 1, 0.3),
  z2 = rnorm(N)
)
df[, x := 0.2 * z1 + 0.3 * z2 + e * 1.5][
    ,y := x + e]

# %%
m_felm = felm(y ~ 1 | 0 | (x ~ z1 + z2), df)
# KP F stat in LFE object / identical in stata's ivreg2
m_felm$stage1$rob.iv1fstat

# %%
m_feols = feols(y ~ 1 | x ~ z1 + z2, df)
fitstat(m_feols, "kpr")
# error

Effective F statistic

Also I was wondering if there's any interest in including the effective F statistic (which is identical to K-P in just-identified single-instrument+single-treatment settings - see [2] for details ).

As part of a review paper I'm working on with coauthors, I'm implementing the Pflueger-Montiel Olea Effective F statistic [1], which is recommended by the review paper by Andrews, Stock, and Sun [2] manually for now using the simplified formula

from the Andrews et al paper (the original paper is written in archaic notation with y and Y as the outcome and treatment respectively; ew).
Pretty straightforward to implement / I can send a pull request with the implementation if you think there's enough interest.

happy to move to related F-stat thread if preferred #117

[1] http://www.joseluismontielolea.com/Montiel-OleaJBES.pdf
[2] https://scholar.harvard.edu/files/wirev_092218-_corrected_0.pdf

[lrberge]

Hey, thanks for the report!

KP problem

I've corrected the bug which was due to the fit method not handling variables with non-standard names.

Now the KP works:

fitstat(m_feols, "kpr")
#> Kleibergen-Paap: stat = 0.673883, p = 0.713951, on 2 DoF.

I must say that the current KP implementation is the same as Stata only when n_inst = n_exo and it diverges when n_exo > n_inst. I don't know the cause and I think I may just put NA instead when n_inst < n_exo. But I would be very happy if you could (or if you know someone who could) fix the KP code (which is a nightmare btw! :-)).

F stat

Isn't it similar to the Wald in fitstat? If you pick the right VCOV, you obtain the same tests as a KP with a single instrument:

m_felm$stage1$rob.iv1fstat
#> $x
#>            p         chi2          df1          p.F            F          df2 
#> 1.452699e-09 4.069969e+01 2.000000e+00 2.176714e-09 2.034984e+01 9.970000e+02 
#> attr(,"formula")
#> ~z1 | z2
#> <environment: 0x00000254374e6460>

fitstat(m_feols, "ivwald", se = "heter")
#> Wald (1st stage), x: stat = 20.3, p = 2.177e-9, on 2 and 997 DoF, VCOV: Heteroskedasticity-robust.

[etiennebacher]

Hello,

I think this bug is back but with another error message. With the Github version:

library(data.table)
library(fixest)

set.seed(42)
N = 1e3
df = data.table(
  e = rnorm(N),
  z1 = rbinom(N, 1, 0.3),
  z2 = rnorm(N)
)
df[, x := 0.2 * z1 + 0.3 * z2 + e * 1.5][
  ,y := x + e]

m_feols = feols(y ~ 1 | x ~ z1 + z2, df)
fitstat(m_feols, "kpr")

Error in kronv %*% vhat : non-conformable arguments

[lrberge]

Indeed, thanks: that's fixed.

[adamaltmejd]

fitstat(x, "kpr") does not seem to work with an estimation run with lean = TRUE even if I manually set vcov to use the same clusters as was initially used by feols(). Is this on purpose? Does fitstat(x, "kpr") need to recalculate se's?

[lrberge]

You need the VCOV indeed.

[adamaltmejd]

Thanks!

[adamaltmejd]

Found another bug. When there are missing values in the instrumented variable the KPR calculation throws:
Error in svd(theta) : infinite or missing values in 'x'

[adamaltmejd]

@lrberge did you look at this?

[adamaltmejd]

I also get different errors when I try to extract KPR stats where the instrument is interacted. Eg:

library(data.table)
library(fixest)
set.seed(42)
N = 1e3
df = data.table(
  e = rnorm(N),
  z1 = rbinom(N, 1, 0.3),
  z2 = rnorm(N),
  int = c(rep(1, 1e3), rep(0, 2e3))
)
df[, x := 0.2 * z1 + 0.3 * z2 + e * 1.5][
  ,y := x + e]

m_feols = feols(y ~ 1 | x * int ~ z1 * int, df)
fitstat(m_feols, "kpr")

Error in solve.default(t(Gmat)) : 
  system is computationally singular: reciprocal condition number = 1.26152e-16

In my actual data I also get the error:

Error in chol.default(crossprod(X_proj)) :                                                                                                                                               
  the leading minor of order 1 is not positive definite  

[lrberge]

Thanks for reporting:

• I won't fix the bug since I don't understand the underlying method well.
• I'm really considering removing this statistics altogether since it's buggy and I can't fix.

On my side I already spent a lot of time implementing it and can't spend more time on it especially that I don't use it. There is just too many implicit knowledge in the paper and I don't have the time to dig into it.
Maybe one of you could take a stab at it? @apoorvalal? you seem super algebra-aware? ;-)

[adamaltmejd]

Makes sense. For the one single instrument setting it does seem to work though, so an option could be to keep that? And throw an error otherwise saying that other specifications are not implemented?

At least when I've been presenting results it is the test of weak instruments that most people have been asking for. I don't think I would be capable of implementing it though...

[adamaltmejd]

I guess it would actually make more sense to implement the Montiel Olea-Pflueger (2013) (weakivtest in Stata) instead as @apoorvalal is suggesting. Then, users could just be referred to use that test instead of KPR.

Did you ever start working on an implementation @apoorvalal? Would be happy to help out if there is anything I can do?

[apoorvalal]

I have a very minimal implementation of it in the package for our paper https://github.com/apoorvalal/ivDiag/blob/master/R/effF.R
Although I've come around to the position that we should probably just be reporting Anderson-Rubin intervals anyway, which is fairly standard in its implementation https://github.com/apoorvalal/ivDiag/blob/master/R/AR_test.R and shouldn't be too bad to incorporate into fixest
(I changed the original implementation from ivpack to work with felm objects [I started work on this before I started using fixest].

[lrberge]

@adamaltmejd: sure I can make it available only for the single case, that seems reasonable.

Otherwise, I'm overhauling completely how fitstats are computed internally. After the changes, making contributions to fit statistics will be super super easy (I'll also write a contributing guide). So if you feel to contribute, it will be easy in a few weeks. Another feature that will be automatic is that when lean = TRUE, all fit stats available for the model will be directly integrated, making other functions work properly down the road.

[adamaltmejd]

That sounds like a great idea. Give me a shout out when its done and I might have time to take a look at implementing F^eff.