Parameter "subset" throws errors in feols() when used with clustered SEs

[noahmbuckley]

Hello, thank you for the excellent package! Very useful and powerful. I've been running some boring, typical feols() regressions with two-way fixed effects and associated clustered SEs. I'm used to being able to subset my data in the regression itself when I use lm() or lmer() or other common functions. When I try that here, I get a variety of errors. Currently I mostly get this error:

Error in vcov.fixest(object, se = se, cluster = cluster, dof = dof, forceCovariance = forceCovariance, :
The length of the clusters (310017) does not match the number of observations in the estimation (149218).

The function call I'm using is of this sort:
mod <- feols(outcome ~ x1 + x2 | fe1 + fe2, data=dat, cluster=c("fe1", "fe2"), subset=(x3==4))

I could send you some data and code if that would help. It happens whether I define a vector (vec <- x3==4) and then include that with subset=vec or whether I do the subsetting right in the feols() call. I think that before I was getting different errors too when I was subsetting directly in the function, although I can't reproduce that right now for some reason.

Thank you!!

[lrberge]

Hi, thanks for reporting.

Unfortunately, I can't replicate the behavior:

base = iris
names(base) = c("y", "x1", "x2", "x3", "fe1")
base$fe2 = sample(1:5, 150, TRUE)

est = feols(y ~ x1 + x2 | fe1 + fe2, base, cluster = c("fe1", "fe2"), subset = ~x3>1)
est
#> OLS estimation, Dep. Var.: y
#> Observations: 93 
#> Fixed-effects: fe1: 2,  fe2: 5
#> Standard-errors: Two-way (fe1 & fe2) 
#>    Estimate Std. Error t value Pr(>|t|)    
#> x1 0.307170   0.068668  4.4733 0.140013    
#> x2 0.849757   0.080917 10.5020 0.060439 .  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> RMSE: 0.318971     Adj. R2: 0.722443
#>                  Within R2: 0.651523

Could you provide a reproducible example?

[noahmbuckley]

Hi Laurent, thanks very much for replying and taking a look.

I found a few things going on. Reproducible code below.

First, basically, it turned out I was passing x3 as an Nx1 matrix instead
of a vector -- this was causing the error. lm() and other functions coerce
that matrix into a vector. Granted, your help text does say to pass a
vector, but this is a totally silent and hard-to-diagnose issue, since data
frames can contain variables as matrices without the user ever knowing.
Obviously now that I've identified the problem I can fix my data, but I
think it would be very nice to coerce just like lm() does.

base = iris
names(base) = c("y", "x1", "x2", "x3", "fe1")
base$fe2 = sample(1:5, 150, TRUE)

base <- rbind(base, data.frame("y"=4.6, "x1"=3.0, "x2"=NA, "x3"=0.6,
"fe1"="setosa", "fe2"=5))
base$x2 <- as.matrix(base$x2)
base$x3 <- as.matrix(base$x3)

est = feols(y ~ x1 + x2 | fe1 + fe2, base, cluster = c("fe1", "fe2"),
subset = ~x3>1)  # error
est = feols(y ~ x1 + x2 | fe1 + fe2, base, cluster = c("fe1", "fe2"),
subset = x3>1)  # error

est = feols(y ~ x1 + x2 | fe1 + fe2, base, cluster = c("fe1", "fe2")) # no
error here or in lm()

mod = lm(y ~ x1 + x2, base, subset=x3>1)

The second thing I discovered is that there cannot be any NAs in x3, which
was certainly unexpected behavior for me. This is not the case for other
variables in the regression.

base = iris
names(base) = c("y", "x1", "x2", "x3", "fe1")
base$fe2 = sample(1:5, 150, TRUE)

base <- rbind(base, data.frame("y"=4.6, "x1"=3.0, "x2"=NA, "x3"=NA,
"fe1"="setosa", "fe2"=5))

est = feols(y ~ x1 + x2 | fe1 + fe2, base, cluster = c("fe1", "fe2"),
subset = ~x3>1)  # error

est = feols(y ~ x1 + x2 | fe1 + fe2, base, cluster = c("fe1", "fe2"))

[lrberge]

Thanks for coming back to me! It turns out that this subset argument was rife with bugs.

I'm coming back to you soon.

[noahmbuckley]

Thank you for being responsive and helpful, good to hear that we can get some bugs squashed. Let me know if I can help on testing or anything.

…

On Tue, Apr 20, 2021 at 6:12 PM Laurent Bergé ***@***.***> wrote: Thanks for coming back to me! It turns out that this subset argument was rife with bugs. I'm coming back to you soon. — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#129 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADSYBY5KI2HV5PCIM3SHCVDTJX36HANCNFSM43AI3SGA> .

[lrberge]

The new version should fix the bug and the request.
I'm not a big fan of automatic conversions of types (since it can be the sign that something went wrong in data creation beforehand), but I've done it for subset, weights and offset, so they now will also accept matrices. Just for you ;-)

[lrberge]

I tend to be very very conservative (code-wise!). But maybe I'm wrong in my approach... I may extend these conversions to more arguments and types (like single column data.frames).

Thanks for the suggestion! And thanks for the support and bug filing!

[noahmbuckley]

Thanks very much, these changes will be extremely helpful in the project I'm currently working on! Again, I appreciate your responsiveness and your providing such a great package. Yes, I can certainly see your point there. I'm no expert on coding best practices, but my sense is that having behavior be a) expected and b) matching commonly used other regression functions like lm() or lmer() or plm() far and away outweighs your preferences on what sort of data people feed to your code. After all, it's the analyst's job to provide good data, not your job to police what my data should look like :) The benefit you get from allowing people to run your code alongside lm() etc without any modification seems to me to be quite large. My view is that it would be best for your function, parameter definitions, etc to match those hyper-popular other functions as closely as possible. That seems like the most conservative coding approach of all. My two cents! Far be it from me to demand anything from someone providing such a great free service, but I would suggest that you include as informative as possible error messages when your code *does* diverge from commonly-expected behavior. Then at least the user would have an idea where to look when your code rejects something like a matrix. Thanks again!

…

On Fri, Apr 23, 2021 at 3:26 AM Laurent Bergé ***@***.***> wrote: I tend to be very very conservative (code-wise!). But maybe I'm wrong in my approach... I may extend these conversions to more arguments and types (like single column data.frames). Thanks for the suggestion! And thanks for the support and bug filing! — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#129 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADSYBY77MY6OGIR5HHONNIDTKEOI5ANCNFSM43AI3SGA> .

[lrberge]

You raise a very interesting point. Why do I keep developing the package, and for whom? I actually don't have a clear answer to it.
Regarding the beginning of the package, the answer is easy: I started its development because I needed it for my research, then I shared it through a package thinking other persons may find it useful.
Regarding its current state it's more blurry: is it for the users or for my research purposes? Although I did add things to the package that I use in my research, many recent additions concern things that I don't (most prominent example is IV estimations). So it seems that I develop the package for the users but I'm still part of the equation.

having behavior be a) expected and b) matching commonly used other regression functions like lm() or lmer() or plm() far and away outweighs your preferences on what sort of data people feed to your code.

Please remember that I still use my software. I don't try to police any one, I simply implement capabilities that fit my workflow, and I'm kind of a control freak. Thus by adding the type conversions a la lm() I end up strictly worse off in my workflow because I would like to be informed if I have non-vectors in my data set. Sure I'm only one person, but why should I implement things that make me worse off?

I would suggest that you include as informative as possible error messages when your code does diverge from commonly-expected behavior.

This comment makes my heart bleed, given the enormous amount of effort that I've put into error handling. Here are the kind of error messages that you can have when using the package:

base = iris ; names(base = c("y", "x1", "x2", "x3", "species")
base$s = as.matrix(base$s)

feols(y ~ x1, base, split = base$s)
#> Error in feols(y ~ x1, base, split = base$s) :
#>   If not a formula, argument 'split' must be either: i) a vector, or ii) a character scalar. Problem: it is not
#> a vector nor a scalar (instead it's a matrix).

base$sub_df = as.data.frame(base$species)
feols(y ~ x1, base, subset = base$sub_df)
#> Error in feols(y ~ x1, base, subset = base$sub_df) : 
#>   If not a formula, argument 'subset' must be either: i) an integer matrix (or vector), or ii) a logical
#> vector. Problem: it is not a matrix nor a vector nor a vector (instead it is a data.frame).

It looks pretty clear to me... The fact that you didn't get such a message means that it was... a bug. Please keep in mind that we're talking about tens of thousands of lines of code, so bugs will always exist!

I suggest you to mess up with feols arguments and look at the error messages that you get; and then do the same exercise with lm/lmer. I'm not sure the comment will still hold ;-).

Just for the record here's an equivalent lmer error message:

fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy, subset = as.data.frame(sleepstudy$Reaction))
#> Error in xj[i] : invalid subscript type 'list'

---

In any case, thanks a lot for filing the issue! The software got better!

[noahmbuckley]

Yes, it does seem like a tricky, fuzzy situation to be in! I see people on academic Twitter talking about the package all the time and, as I'm sure you know, it has become quite popular! Which is great, and again, you're doing a wonderful service. But it is tricky to see how you fit in and who the package is for. I don't know how best to approach the interesting issues and dilemmas you raise. Do you try to find a developer to take over maintenance and development of the package for public consumption? (And you could fork it for yourself to keep things how you like them.) I don't know how possible that is, but that would be one approach, I suppose. I know some package maintainers have done that -- asked for help from the wider community. Very sorry - I didn't mean to imply that you've done a bad job with error handling or anything. It certainly seems like a massive job to keep all of this code working! All the best and thanks for a nice conversation!

…

On Fri, Apr 23, 2021 at 4:40 PM Laurent Bergé ***@***.***> wrote: You raise a very interesting point. Why do I keep developing the package, and for whom? I actually don't have a clear answer to it. Regarding the beginning of the package, the answer is easy: I started its development because I needed it for my research, then I shared it through a package thinking other persons may find it useful. Regarding its current state it's more blurry: is it for the users or for my research purposes? Although I did add things to the package that I use in my research, many recent additions concern things that I don't (most prominent example is IV estimations). So it seems that I develop the package for the users but I'm still part of the equation. having behavior be a) expected and b) matching commonly used other regression functions like lm() or lmer() or plm() far and away outweighs your preferences on what sort of data people feed to your code. Please remember that I still use my software. I don't try to police any one, I simply implement capabilities that fit my workflow, and I'm kind of a control freak. Thus by adding the type conversions *a la* lm() I end up strictly worse off in my workflow because *I* would like to be informed if I have non-vectors in my data set. Sure I'm only one person, but why should I implement things that make me worse off? I would suggest that you include as informative as possible error messages when your code *does* diverge from commonly-expected behavior. This comment makes my heart bleed, given the enormous amount of effort that I've put into error handling. Here are the kind of error messages that you can have when using the package: base = iris ; names(base = c("y", "x1", "x2", "x3", "species")base$s = as.matrix(base$s) feols(y ~ x1, base, split = base$s)#> Error in feols(y ~ x1, base, split = base$s) :#> If not a formula, argument 'split' must be either: i) a vector, or ii) a character scalar. Problem: it is not#> a vector nor a scalar (instead it's a matrix). base$sub_df = as.data.frame(base$species) feols(y ~ x1, base, subset = base$sub_df)#> Error in feols(y ~ x1, base, subset = base$sub_df) : #> If not a formula, argument 'subset' must be either: i) an integer matrix (or vector), or ii) a logical#> vector. Problem: it is not a matrix nor a vector nor a vector (instead it is a data.frame). It looks pretty clear to me... The fact that you didn't get such a message means that it was... a bug. Please keep in mind that we're talking about tens of thousands of lines of code, so bugs will always exist! I suggest you to mess up with feols arguments and look at the error messages that you get; and then do the same exercise with lm/lmer. I'm not sure the comment will still hold ;-). Just for the record here's an equivalent lmer error message: fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy, subset = as.data.frame(sleepstudy$Reaction))#> Error in xj[i] : invalid subscript type 'list' ------------------------------ In any case, thanks a lot for filing the issue! The software got better! — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#129 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADSYBYZL2FT2OP5TCFMPY7DTKHLMDANCNFSM43AI3SGA> .