FEATURE REQUEST: vertexLmer like mincLmer

[gdevenyi]

mincLmer is available for linear mixed effect models, it would be nice of there were a vertex version as well.

[gdevenyi]

As promised, the code we've been using:

### Function for running vertex-wise linear mixed effects models, using lmerTest package. 
### FDR correction included ("step-up" p values in p.adjust).
### Change REML=FALSE to use ML instead of REML, if something starts crapping out for some reason.
### Sample command: left_thickness_results <- vertexLmer(gf, 'y ~ Age + Sex + group*Random_1 + group*Random_2 + (1|Random_1) + (1|Random_2)', left_thickness)
### Matt Park

detach_package <- function(pkg, character.only = FALSE)
{
  if(!character.only)
  {
    pkg <- deparse(substitute(pkg))
  }
  search_item <- paste("package", pkg, sep = ":")
  while(search_item %in% search())
  {
    detach(search_item, unload = TRUE, character.only = TRUE)
  }
}

detach_package(lme)
detach_package(lmerTest)
detach_package(pbkrtest)
detach_package(lme4)

library(lmerTest)

###Degrees of freedom and p values calculated based on Scatterthwaite's approximations

vertexLmer <- function(glim.matrix, statistics.model, vertex.table) {
    number.vertices <- nrow(vertex.table)

    y <- vertex.table[1, ]
    s <- coef(summary(lmer(formula(statistics.model), REML=FALSE, data=glim.matrix)))

    number.terms <- nrow(s)
    variable.names <- rownames(s)

    value <- matrix(data = 0, nrow = number.vertices, ncol = number.terms)
    std.error <- matrix(data = 0, nrow = number.vertices, ncol = number.terms)
    df.value <- matrix(data = 0, nrow = number.vertices, ncol = number.terms)
    t.value <- matrix(data = 0, nrow = number.vertices, ncol = number.terms)
    p.value <- matrix(data = 0, nrow = number.vertices, ncol = number.terms)
    q.value <- matrix(data = 0, nrow = number.vertices, ncol = number.terms)
    r.value <- matrix(data = 0, nrow = number.vertices, ncol = number.terms)

    f<- formula(statistics.model)
    modulo <- 500

    cat("Percent done: ")
    for (v in 1:number.vertices) {
        y <- vertex.table[v, ]
        s= coef(summary(lmer(f, REML=TRUE, data=glim.matrix)))

        value[v, ] <- s[, 1]
        std.error[v, ] <- s[, 2]
        df.value[v, ] <- s[, 3]
        t.value[v, ] <- s[, 4]
        p.value[v, ] <- s[, 5]

        if (v%%modulo == 0) {
            cat(format((v/number.vertices) * 100, digits = 3))
            cat("%  ")
        }
    }

    cat("\n")

    for (i in 1:number.terms) {
        q.value[,i] <- p.adjust(p.value[,i], "fdr")
    }

        for (i in 1:number.terms) {
        r.value[,i] <- (t.value[,i]) / sqrt(((t.value[,i])^2) + df.value[,i])
        }

    colnames(value) <- variable.names
    colnames(std.error) <- variable.names
    colnames(df.value) <- variable.names
    colnames(t.value) <- variable.names
    colnames(p.value) <- variable.names
    colnames(q.value) <- variable.names
    colnames(r.value) <- variable.names

    results <- list(value = value, std.error = std.error, df.value=df.value, t.value = t.value, p.value= p.value, q.value=q.value, r.value=r.value)
    return(results)
}

vertexLmerFDR <- function(results) {
    output_matrix<- data.frame()

    thresholds <- c(0.01,0.05,0.10,0.15,0.20)
    for (i in thresholds){
        threshold=i
        t_results<- matrix(data=NA, nrow=1, ncol=(ncol(results$q.value)-1))

        for (j in 2:ncol(results$q.value)){
            columns<- cbind(results$q.value[,j], results$t.value[,j])
            thresholded<- subset(columns, columns[,1] < threshold)
            t_threshold<- round(suppressWarnings(min(abs(thresholded[,2]))), digits=4)
            if (nrow(thresholded)==0){
                t_threshold<- "NA"
            }
            t_results[1,(j-1)]<- as.character(t_threshold)
        }

        output_matrix<- rbind(output_matrix, t_results)
    }
    colnames(output_matrix)<- colnames(results$q.value)[2:ncol(results$q.value)]
    rownames(output_matrix)<- thresholds
    print(output_matrix)    
}

[cfhammill]

Finished writing our version with parallelization, will be in develop in the next couple house, hoping for a release within the next week or so.