Hi,
given that my first attempts to implement panel logit with so-called
random effects have been 95% but not 100% successful I have now written
a wrapper function to outsource this task to R. Given that this is a
potentially computational-intensive task I think it makes sense to incur
the overhead of invoking R.
If somebody likes to test and feed back the results, they're welcome
(see attached).
You need to have R installed and within R you need the "lme4" package to
be present.
The function interface from a gretl user side is pretty straightforward,
you just specify the dependent variable and then the right-hand side in
a list "X". If you want to let lme4's glmer() do Gauss-Hermitian
quadrature you need to specify the number of points in "quadpoints" (or
you set it to zero to let me choose the default number 32), otherwise
with the default 1 it uses a "Laplacian" quadrature:
RElogitR(series y, const list X, int quadpoints[0::1])
The function returns a bundle with all kind of stuff, but this is far
from finished.
Note that in the output you will see from the R side, the jargon from
the lme4 package calls some things "fixed effects", which is probably a
pretty accurate description, but that we in econometrics usually
wouldn't call fixed effects. It has nothing to do with (conditional)
Fixed-Effects-Logit from gretl's fe_logit package for example. It is
really just the standard explanatory variables without the "random effects".
So let's see if this thing works for you.
cheers,
sven
# Random-Effects Logit outsourced to R
# (glmer function in lme4 package)
# (test code at the bottom)
# Sven Schreiber
function bundle RElogitR(series y, const list X, \
int quadpoints[0::1])
# interpretation of quadpoints:
# default 1 invokes the glmer default of Laplacian,
# setting it to zero invokes our default of 32 for GHQ,
# and any other positive number is passed verbatim.
# a constant term is not explicitly added!
# This function needs R, and the "lme4" package/library
# must be installed inside R.
# set echo off
# set messages off
if !isdummy(y)
funcerr "Dep. var must be a dummy!"
endif
catch series groups = $unit
if $error
funcerr "Is this really a panel dataset?"
endif
# pass the quadpoints option to R
quadpoints = quadpoints==0 ? 32 : quadpoints
matrix mquadp = {quadpoints}
mwrite(mquadp, "RElogitR_mquadp.mat", 1)
# prepare the data for passing to R
list pass = y X groups
foreign language=R
q <- gretl.loadmat("RElogitR_mquadp.mat")
library(lme4)
# use all vars in the model except "groups" (and "y"),
# and include random effects based on "groups"
m <- glmer(y ~ . - groups + (1|groups), data=gretldata,
family="binomial", nAGQ=q)
summary(m)
# retrieve output/results
rbeta <- as.matrix(getME(m, "beta"))
rsigma <- as.matrix(getME(m, "sigma"))
rnumofRE <- as.matrix(getME(m, "n_rtrms"))
# (still missing is VCOV, see merMod ...)
gretl.export(rbeta)
gretl.export(rsigma)
gretl.export(rnumofRE)
end foreign --send-data=pass
bundle b
b.yname = argname(y)
b.Xnames = varname(X)
b.quadpoints = quadpoints
b.beta = mread("rbeta.mat", 1)
b.sigma = mread("rsigma.mat", 1)
b.numofRE = mread("rnumofRE", 1)
return b
end function
## test code
test = 1
if test
open abdata.gdt # --preserve
series binary = (INDOUTPT > INDOUTPT(-1))
list rhs = const EMP WAGE
bundle myRElogit = RElogitR(binary, rhs)
endif
