conducting multi-level regression on ordinal DVs with imputed data in R Do you know of an approach/package that facilitates mixed model regression of ordinal dependent variables on multiply imputed datasets in R?
Ideally, the function takes:
a list of multiply imputed datasets
a list of target variables (dependent variables, one or more)
a list of factors (independent variables, one or more)
a list of dummy coded conditions (for analysis of multifactor IVs)
and returns a table similar to the result of other regressions in r
After extensive searching, I had to create such a function using CLMM in the ordinal package. If you can't answer the first question perhaps you can advise me wrt code adaptations, statistical appropriateness of my approach, efficiency (it takes a LONG time with many imputed datasets), etc
here's some data that mirrors mine
numimpdatasets = 3; N = 170; datalist = list()
for (datasetnum in 1:numimpdatasets){
    dvone = sample(1:5, N, replace=T)   
    dvtwo = sample(1:5, N, replace=T)
    teacher = c('Tom','Dick','Harry')[sample(1:3, N, replace=T)]
    studentclass = c('Class1','Class2','Class3','Class4','Class5',
                     'Class6')[sample(1:6, N, replace=T)]
    aptitude = runif(N, -3.5, 3.5)
    randord = sample(1:3, N, replace=T)
    conddummycode1 = c('a_cond1', 'b_cond2', 'c_cond3')[randord]
    conddummycode2 = c('c_cond1', 'b_cond2', 'a_cond3')[randord]
    datalist[[datasetnum]] = data.frame(cbind(dvone, dvtwo,teacher,
            studentclass, aptitude, conddummycode1, conddummycode2))
    }
dvs = colnames(datalist[[1]])[1:2]
conditions = colnames(datalist[[1]])[6:7]
ivs = c("+as.numeric(aptitude)","+(1|teacher/studentclass)")

 A: here is my approach. I use the coefficient/stderror pooling from the mi package and then calculate the p-value for each factor
Right now I only have the summary table displaying the factors in the regression. The ideal approach would output a table identical to that of summary(clmmModel) but I haven't done the research/found the math for pooling the other output components
require(ordinal)
require(mi)

iterclmm <- function(dvs, ivs, datalist, conditions){       
    for (dv in dvs){
        # # create name for list of target variable models
        varname = paste0(dv,"_modlist")
        ivlist = ""
        for (iv in ivs){ivlist = paste0(ivlist,iv)}
        # # roll through each dummy-coded condition designation
        for (condition in conditions){
            templist = list()
            # # this count used to create a list of models/dataset
            count = 1
            # # iterate through each dataset and run model
            for (dataset in datalist){
                # # set up  clmm regression
                tempcom = paste0('model <- clmm(as.factor(',dv,') ~ as.factor(',condition,')',ivlist,',data = dataset)')
                # # run it
                eval(parse(text=tempcom))
                # # save it, print model summary
                eval(parse(text=paste0("templist[[",count,"]] = model")))
                count = count + 1
                }   
        # # # # create pooled coefficient table
            modelrownames = rownames(coefficients(summary(templist[[1]])))
            coeflist = list()
            # # HACK TO MAKE LIST WORK - NEED BETTER SOLUTION FOR ADDING TO EMPTY LIST
            coeflist[[1]] = 0
            # # make a list of coefficients from all models
            for (model in templist){
                coeflist[[length(coeflist)+1]] = coefficients(summary(model))           
                }
            coeflist = coeflist[2:length(coeflist)]
            # # transform coefficient list
            coefest = list(); stderrors = list() 
            zvals = list(); pvals = list()
            count = 1
            # # assign vals 
            for (coeftable in coeflist){
                temptable = t(coeftable)
                coefest[[count]] = temptable[1,]; stderrors[[count]] = temptable[2,]
                zvals[[count]] = temptable[3,]; pvals[[count]] = temptable[4,]
                count = count + 1
                }
            tempvals = data.frame(cbind(coefest, stderrors, zvals, pvals))
            eval(parse(text=paste0(dv, "_regvals <<- tempvals")))
            # # create regression results table for target variables
            print(dv)
            print(condition)
            coefsandstderror = mi.pooled(coefest, stderrors)
            coef.est = coefsandstderror$coefficients; std.error = coefsandstderror$se
            zvalues = coef.est / std.error; pvalues = 2 * (1- pnorm(abs(zvalues)))
            pooltable = data.frame(cbind(coef.est,std.error,zvalues,pvalues))
            eval(parse(text=paste0(dv, "_sum <<- pooltable")))
            print(pooltable)
            cat(rep("\n",2))
            }       
        cat(rep("\n",4))
        }
    }

