Edit: I have updated this question with the solution below:

I am trying to understand how pglm estimates fixed effects models for binomial dependent variables by replicating the same model with the dummy variable method and the variance decomposition method in R. I am able to replicate the linear models produced via plm perfectly, but I cannot figure out how to replicate pglm's estimates via the same procedures. Would anyone know what I am missing or doing wrong here?

Below illustrates how I replicate plm and pglm using linear models and nonlinear models with plm's "Produc" dataset. If the required libraries are installed, The script should work if copied and pasted into directly into one's R console.

#0) set seed

#1) load libraries

#2) prepare data
    #a) prepare base data
        #i) load produc
            data("Produc", package = "plm")
        #ii) convert to data table
            Produc <- data.table(Produc)
        #iii) create some fake dummy variables, rescale independent variables, remove the rest of the df
            df <- Produc[,.(
              binvar = rbinom(nrow(Produc),1,.02),
              dv.binvar = rbinom(nrow(Produc),1,.04 * abs(scale(pcap)) * abs(scale(emp))),
              pcap = pcap/10000,
              pc = pc/10000,
              emp = emp/10000,
              state = state
        #iv) make binvar time-invariant
            df[, ("binvar") := lapply(.SD, max), by = "state", .SDcols = "binvar"]

    #b) prepare variance decomposition version of data
        #i) find mean difference for all variables
            df.decomp <- df[, lapply(.SD, function(x) x - mean((x))), .SDcols = c("pcap", "pc", "emp", "binvar", "dv.binvar"), by = "state"]
        #ii) add non-transformed versions of dummy variables back into the data
              ti.binvar = df$binvar,
              ti.dv.binvar = df$dv.binvar

#3) linear models
    #a) estimate linear models
        m.plm <- plm(dv.binvar ~ pcap + pc + emp + emp:binvar, data = df, index = c("state"))
        m.lm <- lm(dv.binvar ~ pcap + pc + emp + emp:binvar + state - 1, data = df)
        m.decomp <- lm(dv.binvar ~ pcap + pc + emp + emp:ti.binvar - 1, data = df.decomp)

    #b) format results of linear models 
        r.plm <- round(summary(m.plm)$coefficients,2)
        r.lm <- round(summary(m.lm)$coefficients,2)
        r.lm <- r.lm[!grepl("state", row.names(r.lm)),]
        r.decomp <- round(summary(m.decomp)$coefficients,2)

    #c) print results of linear models
            #all estimates are the same

#4) nonlinear models
    #a) estimate nonlinear models
        #m.pglm <- pglm(dv.binvar ~ pcap + pc + emp + emp:binvar, data = df, index = c("state"), family = binomial)
        m.bife <- bife(dv.binvar ~ pcap + pc + emp + emp:binvar | state, data = df, model = "logit")
        m.glm <- glm(dv.binvar ~ pcap + pc + emp + emp:binvar + state - 1, data = df, family = binomial)
        m.gdecomp <- glm(ti.dv.binvar ~ pcap + pc + emp + emp:ti.binvar - 1, data = df.decomp, family = binomial)

    #b) format results of nonlinear models
        r.bife <- round(summary(m.bife)$cm,2)
        r.glm <- round(summary(m.glm)$coefficients,2)
        r.glm <- r.glm[!grepl("state", row.names(r.glm)),]
        r.gdecomp <- round(summary(m.gdecomp)$coefficients,2)

    #c) print results of nonlinear models
            #The variance decomposition method is very different!

Thanks, for all of your feedback. It seems I made two key errors:

  • Computational implementation: I assumed that the default method in pglm is "within", as it is in plm. Moreover, it appears that pglm does not even support within person fixed effects! I have updated the syntax in my question so that it uses bife. This returns the same estimates as my dummy variable method, but not my variance decomposition model.

  • The variance decomposition or "demeaning" method cannot work because the model is nonlinear. Demeaning variables wholly changes their relationship with the dependent variable. Moreover, demeaning the dependent variable expands its range beyond 0-1.

It seems that the best method for implementing logistic fixed effects models in R is to implement a "pseudo-demeaning" algorithm, such as that the one used by bife. A detailed explanation for this method is described in: Fernández-Val, I., & Weidner, M. (2018). Fixed effects estimation of large-t panel data models. Annual Review of Economics, 10, 109-138.

  • $\begingroup$ (Note to any reviewers considering VTC: this is clearly R-heavy, but it does seem like a statistical question to me.) $\endgroup$ Nov 5, 2019 at 19:09
  • $\begingroup$ I see this mainly as a question as to how to implement certain procedures in R and for that reason I do think it is appropriate to close as off topic. $\endgroup$ Nov 5, 2019 at 19:52
  • $\begingroup$ I am mainly concerned with understanding the methods underlying pglm, rather than how to actually code it in R (I can easily figure this out myself if I understand what's going on under the hood). I thought presenting my code would illustrate my understanding of the method, and illustrate how I am wrongly assuming pglm works under the hood, but perhaps this is distracting. Would either of you have any tips for how to more effectively ask my question? $\endgroup$
    – Brian A
    Nov 5, 2019 at 20:50
  • $\begingroup$ if you found my answer useful pls. accept and upvote. $\endgroup$ Nov 7, 2019 at 22:13
  • 1
    $\begingroup$ Thanks Jesper. My stackexchange points are too low to allow my upvote to be reflected in the channel. I have accepted your answer and updated my answer/question above with your correction. $\endgroup$
    – Brian A
    Nov 11, 2019 at 16:33

1 Answer 1


I think you are making two mistakes (assuming I have understood the code correctly). The first mistake is related to software use, the second is perhaps more related to understanding of the statistical model. I start with pointing out the second mistake: You cannot remove fixed effects in a non-linear model simply by demeaning the variables and hence the call

m.gdecomp <- glm(ti.dv.binvar ~ pcap + pc + emp + emp:ti.binvar - 1, data = df.decomp, family = binomial)

fitting a logit model to demeaned variables should not be expected to somehow be an implementation of a logit-panel data model with state-fixed effects (individual fixed effects).

The next mistake is that you are missing the fact that the pglm does not estimate fixed effects models (at least not with the binomial link). If you still want to use it for estimation of models that include fixed effects then you should insert appropriate dummy variable (in you case state)

 m.pglm.dummy <- pglm(dv.binvar ~ pcap + pc + emp + emp:binvar + state - 1, data = df, index = c("state"), family = binomial)

which will amount to the same as

 m.glm <- glm(dv.binvar ~ pcap + pc + emp + emp:binvar + state - 1, data = df, family = binomial)

Why did you makes this errorouness assumption? Well I can only speculate, however when you use plm, then without specifying model as "within" and simply providing the index="state" make plm estimate a one way fixed effect model instead of simply a pooled OLS panel data model. If however you specify model="within" in pglm then you get an error (the estimator is not implemented see pglm fail for within). And if you do not specify "within" then it defaults to model="pooled".

So what should you do? I would take a look at the implementation in the package "bife"

bife can be used to fit fixed effects binary choice models (logit and probit) based on an unconditional maximum likelihood approach. It is tailored for the fast estimation of binary choice models with potentially many individual fixed effects. The routine is based on a special pseudo demeaning algorithm derived by Stammann, Heiss, and McFadden (2016). The estimates obtained are identical to the ones of glm, but the computation time of bife is much lower.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.