UP | HOME

Validity of using random-effects parameter estimates from nonconverged models

The question here was whether the random effect parameter estimates from a nonconverged model in lmer are informative in deciding which higher-order term to remove. This analysis considered only within-item designs, where there were two random slopes that could be ranked in terms of their relative size. The correlation between the rank of the estimates from maximal-random-effects models that did not converge and the rank of the generative parameters for the 12-item design was .74, with \(p<.0001\) (N=387 nonconverging models); for the 24-item design it was .83, with \(p<.0001\) (N=192 nonconverging models). Thus, using the nonconverged estimates seems a good strategy.

Source code

library(simgen)

get.ncvix <- function(x) {
  # return a logical vector indicating which runs did not converge (TRUE)
  (1:nrow(x))[x[,"fm"]>0]
}

param.mx <- createParamMx(h0=FALSE)

# load in the maximal random effects models
load(file="test/fitlmer.h0.wswi.12.RS.RData")
ix.12 <- get.ncvix(mx) # keep only those that didn't converge
x.12 <- mx[ix.12,]
prm.12 <- param.mx[ix.12,]

load(file="test/fitlmer.h0.wswi.24.RS.RData")
ix.24 <- get.ncvix(mx) # keep only those that didn't converge
x.24 <- mx[ix.24,]
prm.24 <- param.mx[ix.24,]

# recreate the dataframes that the models will be fit to
df.12 <- mclapply(1:nrow(prm.12),
                  function(x) {mkDf(nsubj=24,nitem=12,params=prm.12[x,],wsbi=FALSE)})

# refit the models
library(lme4)
lmer.12 <- mclapply(df.12,
                    function(x) {
                      lmer(Resp ~ Cond + (1 + Cond | SubjID) + (1 + Cond | ItemID), data=x, na.action=na.omit, REML=FALSE)
                    })

# recreate the dataframes that the models will be fit to
df.24 <- mclapply(1:nrow(prm.24),
                  function(x) {mkDf(nsubj=24,nitem=24,params=prm.24[x,],wsbi=FALSE)})

# refit the models
lmer.24 <- mclapply(df.24,
                    function(x) {
                      lmer(Resp ~ Cond + (1 + Cond | SubjID) + (1 + Cond | ItemID), data=x, na.action=na.omit, REML=FALSE)
                    })

getlmerrank <- function(x) {
  vc <- VarCorr(x)
  myrank <- ifelse(abs(vc$SubjID[2,2]) <= abs(vc$ItemID[2,2]),0,1)
  return(myrank)
}

rank.prm.12 <- ifelse(abs(prm.12[,"t11"])<=abs(prm.12[,"w11"]),0,1)
rank.lmer.12 <- unlist(lapply(lmer.12, getlmerrank))

rank.prm.24 <- ifelse(abs(prm.24[,"t11"])<=abs(prm.24[,"w11"]),0,1)
rank.lmer.24 <- unlist(lapply(lmer.24, getlmerrank))

ranks.12 <- cbind(nitem=12, param.rank=rank.prm.12, lmer.rank=rank.lmer.12)
ranks.24 <- cbind(nitem=24, param.rank=rank.prm.24, lmer.rank=rank.lmer.24)

ranks <- rbind(ranks.12,ranks.24)

cor(ranks[ranks[,"nitem"]==12,])
cor(ranks[ranks[,"nitem"]==24,])

shufflemx <- function(x) {
  ix.12 <- (1:nrow(x))[x[,"nitem"]==12]
  ix.24 <- (1:nrow(x))[x[,"nitem"]==24]
  shuf.12 <- cbind(nitem=12,param.rank=x[ix.12,"param.rank"],lmer.rank=x[sample(ix.12),"lmer.rank"])
  shuf.24 <- cbind(nitem=24,param.rank=x[ix.24,"param.rank"],lmer.rank=x[sample(ix.24),"lmer.rank"])
  return(rbind(shuf.12,shuf.24))
}

getcor <- function(x) {
  c12 <- cor(x[x[,"nitem"]==12,c("param.rank","lmer.rank")])[1,2]
  c24 <- cor(x[x[,"nitem"]==24,c("param.rank","lmer.rank")])[1,2]
  return(c(c12=c12,c24=c24))
}

library(multicore)
nmc <- 9999
perm.mx <- matrix(ncol=2, nrow=nmc+1)
perm.mx[1,] <- getcor(ranks)
colnames(perm.mx) <- names(getcor(ranks))
mx2 <- mclapply(1:nmc, function(x) {getcor(shufflemx(ranks))})
perm.mx[2:(nmc+1),] <- matrix(unlist(mx2), ncol=2, byrow=T)
nge <- apply(apply(perm.mx, 2, function(x) {abs(x)>=abs(x[1])}), 2, sum)
print(nge/nrow(perm.mx))
save.image(file="ranking-test-results2.RData")

Author: Dale J. Barr, Roger Levy, Christoph Scheepers, and Harry J. Tily (daleb@daleb-pc)

Date: April 23, 2012

Generated by Org version 7.8.06 with Emacs version 23

Validate XHTML 1.0