R/mpredict.R

In mpplot: Plots for the R Package MCMCpack

mpredict = function(model, newdata, probs=NULL) {
  posterior = data.frame(as.matrix(model))
  numparam  = dim(newdata)[2]
  params    = names(newdata)
  pnames    = colnames(posterior)      # names in posterior
  ynewpred  = posterior[,1]            # intercept
  mdata     = as.character(attr(model, "call"))[3] # dataset used
  datatypes = sapply(newdata, class)
  # datatypes[grep("^factor$", datatypes)] <- "character"
  for(i in 1:numparam) {
    if(datatypes[i] %in% c("factor", "character")) { # categorical
      test = paste0(colnames(newdata[i]), as.character(newdata[[i]]))
      ppos = 1 # will be replaced if matched, here to avoid multiplying by nothing
      if(test %in% pnames) { ## XXX replace with precise grep?
        ppos = grep(paste0("^",test,"$"), pnames)
        ismatch = 1
        } else {
        ismatch = 0
        }
      ynewpred = ynewpred + posterior[,ppos] * ismatch
      } else {
      ppos = grep(paste0("^",params[i],"$"), pnames) # position in posterior
      ynewpred = ynewpred + posterior[,ppos] * as.numeric(newdata[[i]])
      }
    }
  if(dim(data.frame(ynewpred))[1] == 0) {stop("Error: Could not calculate predictions. Does newdata match the model?")}
  # code for interaction terms:
    msum = summary(model)$statistics   # model summary
    mvar = rownames(msum)              # variable names in summary
    mvar = mvar[2:(length(mvar)-1)]    # cut intercept and sigma2
    intec = grepl(":", mvar)           # identify interaction terms
    intew = which(intec, mvar)         # which terms are interaction terms?
    intes = mvar[intew]                # actual interaction terms
    numint = length(intes)             # how many interaction terms are there
  if (numint > 0) {                    # interaction terms are present
    for(m in 1:numint) {              # repeat for each interaction terms
      # identify parts
      parts = unlist(strsplit(intes[m], split=":")) # parts of current interaction term
      lparts = length(parts)
      npos = rep(NA, lparts)           # positions of the interaction terms in newdata
      mval = rep(NA, lparts)            # include (matched)
      catvars = names(datatypes)[datatypes%in%c("character", "factor")] # list of cateorical variables
      catlen  = length(catvars)                                         # number of categorical variables
      for (p in 1:lparts) {            # identify position of each part
        pinnew = any(grepl(paste0("^",parts[p],"$"), params)) # part is in newdata
        pinpos = any(grepl(paste0("^",parts[p],"$"), pnames)) # part is in model (posterior)
        if(pinpos & !pinnew){
          # part in posterior but not newdata: categorical
          if(catlen > 0){
            for(v in 1:catlen){
              tescan = unique(get(catvars[v], pos=get(mdata))) # candidates
              teslen = length(tescan)
              for(w in 1:teslen){
                if(any(grepl(paste0("^",paste0(catvars[v], tescan[w]),"$"), pnames))){ # is in model (with cateogory)
                  npos[p] = grep(paste0("^",catvars[v],"$"), params)
                  mval[p] = ifelse(as.character(newdata[catvars[v]]) == tescan[w], 1, 0)
                  }
                }
              }
            }
          }
          else {
          # continuous
          npos[p] = grep(paste0("^",parts[p],"$"), params) # position in newdata
          mval[p] = as.numeric(newdata[npos[p]])
          }
        }
      ynewpred = ynewpred + posterior[,1+numparam+m] * prod(mval) # interaction part
      }
    }
  mea = mean(ynewpred)
  if (is.null(probs)) {
    ret = mea
    } else {
    names(mea) = "Mean"
    ret = c(mea, quantile(ynewpred, probs=probs))
    }
  return(ret)
}