# 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
ix.12 <- get.ncvix(mx) # keep only those that didn't converge
x.12 <- mx[ix.12,]
prm.12 <- param.mx[ix.12,]

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")


Date: April 23, 2012

Generated by Org version 7.8.06 with Emacs version 23

Validate XHTML 1.0