Probability predictions with model averaged Cumulative Link Mixed Models fitted with clmm in ordinal package

Question

I found that the predict function is currently not implemented in cumulative link mixed models fitted using the clmm function in ordinal R package. While predict is implemented for clmm2 in the same package, I chose to apply clmm instead because the later allows for more than one random effects. Further, I also fitted several clmm models and performed model averaging using model.avg function in MuMIn package. Ideally, I want to predict probabilities using the average model. However, while MuMIn supports clmm models, predict will also not work with the average model.

Is there a way to hack the predict function so that the function not only could predict probabilities from a clmm model, but also predict using model averaged coefficients from clmm (i.e. object of class "averaging")? For example:

require(ordinal)
require(MuMIn)

mm1 <- clmm(SURENESS ~ PROD + (1|RESP) + (1|RESP:PROD), data = soup,
        link = "probit", threshold = "equidistant")

## test random effect:
mm2 <- clmm(SURENESS ~ PROD + (1|RESP) + (1|RESP:PROD), data = soup,
        link = "logistic", threshold = "equidistant")

#create a model selection object
mm.sel<-model.sel(mm1,mm2)

##perform a model average
mm.avg<-model.avg(mm.sel)


#create new data and predict
new.data<-soup

##predict with indivindual model
predict(mm1, new.data)

I got the following error message: In UseMethod("predict") : no applicable method for predict applied to an object of class "clmm"

 ##predict with model average
 predict(mm.avg, new.data)

Another error is returned: Error in predict.averaging(mm.avg, new.data) : predict for models 'mm1' and 'mm2' caused errors

Sterigma

Why isn't this question directed to the package authors? This seems extremely likely to be "too broad" in that it would require both theoretic and implementation exertion to do in a principled manner. – Herein

Answer 1

I've been using clmm as well and yes I confirm predict.clmm is NOT (yet?) implemented. I didn't yet check the source code for fake.predict.clmm. It might work. If it doesn't, you're stuck with doing stuff by hand or using predict.clmm2.

Strophanthin

As it’s currently written, your answer is unclear. Please edit to add additional details that will help others understand how this addresses the question asked. You can find more information on how to write good answers in the help center. – Serene

Answer 2

I found a potential solution (pasted below) but have not been able to make work for my data.

Solution here: https://gist.github.com/mainambui/c803aaf857e54a5c9089ea05f91473bc

I think the problem is the number of coefficients I am using but am not experienced enough to figure it out. Hopefully this helps someone out though.

This is the model and newdata that I am using, though it is actually a model averaged version. Same predictors though.

ma10 <- clmm(Location3 ~ Sex * Grass3 + Sex * Forb3 + (1|Tag_ID), data = 
IP_all_dunes)
ma_1 <- model.avg(ma10, ma8, ma5)##top 3 models
new_ma<- data.frame(Sex = c("m","f","m","f","m","f","m","f"),
                Grass3 = c("1","1","1","1","0","0","0","0"),
                Forb3 = c("0","0","1","1","0","0","1","1"))


# Arguments:
#  - model = a clmm model
#  - modelAvg = a clmm model average (object of class averaging)
#  - newdata = a dataframe of new data to apply the model to
# Returns a dataframe of predicted probabilities for each row and response level
fake.predict.clmm <- function(modelAvg, newdata) {
  # Actual prediction function
  pred <- function(eta, theta, cat = 1:(length(theta) + 1), inv.link = plogis) {
    Theta <- c(-1000, theta, 1000)
    sapply(cat, function(j) inv.link(Theta[j + 1] - eta) - inv.link(Theta[j] - 
eta))
  }

  # Multiply each row by the coefficients
  #coefs <- c(model$beta, unlist(model$ST))##turn off if a model average is used
  beta <- modelAvg$coefficients[2,3:12]
  coefs <- c(beta, unlist(modelAvg$ST))

  xbetas <- sweep(newdata, MARGIN=2, coefs, `*`)

  # Make predictions
  Theta<-modelAvg$coefficients[2,1:2]
  #pred.mat <- data.frame(pred(eta=rowSums(xbetas), theta=model$Theta))
  pred.mat <- data.frame(pred(eta=rowSums(xbetas), theta=Theta))
  #colnames(pred.mat) <- levels(model$model[,1])
  a<-attr(modelAvg, "modelList")
  colnames(pred.mat) <- levels(a[[1]]$model[,1])

  pred.mat
}

Episcopalism