I am analyzing some data that came from demographic health surveys.  The
data contain information for individuals within households, that are
located within clusters, that are located within survey years, that are
located within countries.  We are trying to find the best model from a
subset of predictors, and all models must contain the random variable of
household within cluster within year within country.

We are running models on a server with 64GB memory and 6 CPU cores

Data:  Please email me if you are willing to look at this for us and I will
send you the data via file transfer. It's a big file and I don't want to
post it online just now.

First, we tested a simple linear mixed model using the lmer package:

########## Install Packages ##########

library(lme4)
 library(glmulti)

########## Clear all memory/objects ##########
 rm(list=ls())
 ########## Read in Data ##########

mydata = read.csv("kr.and.GIS.cleaned.Residents.only.csv")

###################################################################
########## Try lmer model NO Interactions ##########
###################################################################
ptm <- proc.time()

lmer.model = lmer(stunt.dhs ~
       dis_ed_des+
tc_pa+
avg_clu_tc+
 dist_road+
pden_lscan+
URBAN_RURA+
time.water+
 wealth.index +
        (1|country.code.short/year/cluster/household),
        mydata,
REML = F)
print(lmer.model)

proc.time() - ptm


This works fine, runs in 12 minutes, and gives the following output:


Linear mixed model fit by maximum likelihood ['lmerMod']
Formula: stunt.dhs ~ dis_ed_des + tc_pa + avg_clu_tc + dist_road +
pden_lscan +      URBAN_RURA + time.water + wealth.index + (1 |
country.code.short/year/cluster/household)
   Data: mydata
      AIC       BIC    logLik  deviance
 218267.5  218420.4 -109116.7  218233.5
Random effects:
 Groups                                        Name        Std.Dev.
 household:(cluster:(year:country.code.short)) (Intercept) 0.5584
 cluster:(year:country.code.short)             (Intercept) 0.3727
 year:country.code.short                       (Intercept) 0.0000
 country.code.short                            (Intercept) 0.3028
 Residual                                                  1.3819
Number of obs: 59405, groups:
household:(cluster:(year:country.code.short)), 40703;
cluster:(year:country.code.short), 7436; year:country.code.short, 27;
country.code.short, 21
Fixed Effects:
        (Intercept)           dis_ed_des                tc_pa
 avg_clu_tc
         -1.324e+00            1.348e-07           -7.806e-04
 -5.180e-04
          dist_road           pden_lscan          URBAN_RURAU
 time.water
         -1.435e-06            9.003e-06            6.637e-02
  6.097e-05
 wealth.indexpoorer  wealth.indexpoorest   wealth.indexricher
wealth.indexrichest
         -7.814e-02           -1.794e-01            1.469e-01
  4.390e-01  > > proc.time() - ptm   user  system elapsed
727.842  13.781 741.053


Next, we tried to do model selection using the genetic algorithm of
glmulti.  We eventually want to include more predictors than in this
current model AND we want to include interactions.  This could lead to over
one billion candidate models.  For this reason, we want to eventually use
the genetic algorithm and NOT all possible subsets.  However, as a first
step, we will use this smaller list of predictors and level = 1 so that no
interactions are considered.

First, we create the lmer wrapper and then just run the model to see how
many candidate models we would have with this situation

###################################################################
########## Create lmer wrapper to paste random variable to all models ######
###################################################################

lmer.glmulti <- function (formula, data, random = "", ...) {
lmer(paste(deparse(formula), random), data = data, REML=F, ...)
}


##################################################################
########## Find canditate number of models NO interactions ##########
###################################################################


candidate <- glmulti(stunt.dhs ~
 dis_ed_des*
        tc_pa*
         avg_clu_tc*
dist_road*
pden_lscan*
 URBAN_RURA*
time.water*
wealth.index,
data=mydata,
 level = 1,
fitfunc = lmer.glmulti,
        random = "+(1|country.code.short/year/cluster/household)",
        method = "d")


With this model, the same that we used for the lmer function, we get the
following error.  We also get this error if we try to run the same code but
with method = "g", and confsetsize=5)



Initialization... Error in
factor(household:(cluster:(year:country.code.short))) :    negative
length vectors are not allowed



If, however, we delete the household part of the random factor, the
function works:


candidate <- glmulti(stunt.dhs ~
 dis_ed_des*
        tc_pa*
         avg_clu_tc*
dist_road*
pden_lscan*
 URBAN_RURA*
time.water*
wealth.index,
data=mydata,
 level = 1,
fitfunc = lmer.glmulti,
        random = "+(1|country.code.short/year/cluster)",

        method = "d")


Gives:


Initialization...
TASK: Diagnostic of candidate set.
Sample size: 59405
2 factor(s).
6 covariate(s).
0 f exclusion(s).
0 c exclusion(s).
0 f:f exclusion(s).
0 c:c exclusion(s).
0 f:c exclusion(s).
Size constraints: min =  0 max = -1
Complexity constraints: min =  0 max = -1
Your candidate set contains 256 models.


And running the following code seems to be working (slowly):

###################################################################
########## Run models with NO Interactions using genetic algorithm #####
###################################################################

ptm <- proc.time()

model <- glmulti(stunt.dhs ~
dis_ed_des*
        tc_pa*
        avg_clu_tc*
        dist_road*
pden_lscan*
 URBAN_RURA*
time.water*
wealth.index,
data=mydata,
 level = 1,
fitfunc = lmer.glmulti,
        random = "+(1|country.code.short/year/cluster)",
        crit = bic,
        imm = 0.5,
        confsetsize = 5)

CAN ANYONE EXPLAIN WHY WE GET THE NEGATIVE LENGTH VECTORS ERROR IF WE
INCLUDE HOUSEHOLD NESTED WITHIN THE OTHER RANDOM FACTORS FOR GLMULTI BUT
NOT FOR LMER?  ANY SUGGESTIONS FOR WHAT WE CAN DO TO DIAGNOSE OR FIX THE
PROBLEM?

Thanks!

-- 
Alicia Ellis
Postdoc
Gund Institute for Ecological Economics
University of Vermont
617 Main Street
Burlington, VT  05405
(802) 656-1046
http://www.wcs-heal.org
http://www.uvm.edu/~aellis5
<http://entomology.ucdavis.edu/faculty/scott/aellis/>

	[[alternative HTML version deleted]]


