Replicating results from R pglm (binomial fixed effects model) 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
    set.seed(1)


#1) load libraries
    library(plm)
    library(data.table)
    library(bife)


#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
            df.decomp[,':='(
              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
        r.plm
        r.lm
        r.decomp
            #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
        r.bife
        r.glm
        r.gdecomp
            #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.
 A: 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. 
