# 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