[R-sig-ME] perturb with lme4

[Austin Frank]

Hello!

I'd like to evaluate the degree of colinearity introduced into my model
by two categorical independent variables.  I'd like to use a
perturbation analysis, and so reached for the 'perturb' package.
Unfortunately, the perturb function as written doesn't work with lmer
models.

The funtion directly accesses the "call" slot of its model argument (frm
<- mod$call$formula, for example).  I attempted the obvious fix of
either using an accessor function (frm <- formula(mod)) or using the
right slotting operator (frm <- mod at call$frm).  This naive replacement
attempt didn't get me very far.

Can anyone suggest a way to modify the perturb function as written, some
work-alike code for an lmer model, or an alternative way for assessing
colinearity introduced by categorical variables in an lmer model?

My (half-hearted) attempt at a modified perturb function is below.

Thanks,
/au


## should be called with
# ptb <- my.perturb(mod.lmer,
#     pfac=list(
#         list("fac1",pcnt=95),
#         list("fac2",pcnt=95)))

my.perturb <- function
( mod, pvars = NULL,
 prange = NULL,
 ptrans = NULL,
 pfac = NULL, 
 uniform = FALSE,
 niter = 100) 
{
    cutsp <- function(indx, tbl) {
        findInterval(runif(1), tbl[indx, ], rightmost.closed = TRUE)
    }
    if (is.null(formula(mod))) 
        stop("First argument does not contain a formula")
    stopifnot(is.list(pfac) || is.null(pfac))
    nms <- all.vars(terms(mod))
    stopifnot(all(pvars %in% nms))
    result <- NULL
    ncases <- length(get(nms[1]))
    frm <- deparse(formula(mod), width.cutoff = 500)
    result$formula <- frm
    allb <- coefficients(mod)
    if (length(pvars) > 0) {
        stopifnot(is.vector(get(pvars)))
        stopifnot(length(pvars) == length(prange))
        result$pvars <- pvars
        result$prange <- prange
        if (length(ptrans) > 0) 
            result$ptrans <- ptrans
        b <- make.names(c(nms, pvars), unique = TRUE)
        pvars.1 <- b[(length(nms) + 1):length(b)]
        for (i in 1:length(pvars)) {
            inp <- paste("\\<", pvars[i], "(\\>[^.]|$)", sep = "")
            outp <- paste(pvars.1[i], "\\1", sep = "")
            frm <- gsub(inp, outp, frm)
            ptrans <- gsub(inp, outp, ptrans)
        }
        result$ptrans2 <- ptrans
    }
    if (length(pfac[[1]]) > 0) {
        rcls.tbl <- NULL
        pfac.1 <- NULL
        if (is.list(pfac[[1]])) 
            n <- length(pfac)
        else n <- 1
        for (i in 1:n) {
            if (n == 1) 
                lstnm <- pfac
            else lstnm <- pfac[[i]]
            stopifnot(all(lstnm[[1]] %in% nms))
            b <- make.names(c(nms, lstnm[[1]]), unique = TRUE)
            pfc <- b[(length(nms) + 1):length(b)]
            inp <- paste("\\<", lstnm[[1]], "(\\>[^.]|$)", sep = "")
            outp <- paste(pfc, "\\1", sep = "")
            frm <- gsub(inp, outp, frm)
            rcls <- do.call("reclassify", lstnm)
            rcls.tbl <- c(rcls.tbl, list(rcls))
            pfac.1 <- c(pfac.1, pfc)
        }
        result$reclassify.tables <- rcls.tbl
    }
    result$formula2 <- frm
    mod at call$formula <- as.formula(frm)
    if (uniform) {
        ranexp <- "runif(ncases,-prange[i],prange[i])"
        result$distribution <- "uniform"
    }
    else {
        ranexp <- "rnorm(ncases,0,prange[i])"
        result$distribution <- "normal"
    }
    for (k in 1:niter) {
        if (length(prange) > 0) {
            for (i in 1:length(prange)) {
                assign(pvars.1[i], get(pvars[i]) + eval(parse(text = ranexp))
             }
        }
        for (trans in ptrans) eval(trans)
        if (length(pfac[[1]]) > 0) {
            for (i in 1:length(rcls.tbl)) {
                tbl <- rcls.tbl[[i]]$cum.reclass.prob
                tbl <- cbind(0, tbl)
                assign("tmpvar", as.numeric(get(rcls.tbl[[i]]$variable)))
                assign(pfac.1[i], sapply(tmpvar, cutsp, tbl))
                assign(pfac.1[i], as.factor(get(pfac.1[i])))
            }
        }
        mod2 <- eval(mod at call)
        allb <- rbind(allb, coefficients(mod2))
    }
    rownames(allb) <- NULL
    result$coeff.table <- allb
    class(result) <- "perturb"
    result
}

-- 
Austin Frank
http://aufrank.net
GPG Public Key (D7398C2F): http://aufrank.net/personal.asc