R/addreg.reparameterise.r

In addreg: Additive Regression for Discrete Data

addreg.reparameterise <- function(nn.coefs, terms, data, type = c("cem", "em"), allref, mono, design.ref)
{
  type <- match.arg(type)
  termlabels <- attr(terms,"term.labels")
  design.type <- sapply(allref, attr, "type")
  ref.vector <- as.vector(design.ref, mode = "integer")
  
  npar.o <- npar.n <- 1L + sum(design.type == 1)
  if (sum(design.type == 2) > 0) {
    npar.o <- npar.o + sum(sapply(sapply(termlabels[design.type == 2], function(x) gsub("`", "", x)), 
                                  function(x) nlevels(factor(data[[x]])) - as.numeric(type == "cem")))
    npar.n <- npar.n + sum(sapply(sapply(termlabels[design.type == 2], function(x) gsub("`", "", x)), 
                                  function(x) nlevels(factor(data[[x]])) - 1))
  }
  if (sum(design.type == 3) > 0) {
    npar.o <- npar.o + sum(sapply(sapply(termlabels[design.type == 3], function(x) gsub("`", "", x)), 
                                  function(x) nlevels(factor(data[[x]])) - 1))
    npar.n <- npar.n + sum(sapply(sapply(termlabels[design.type == 3], function(x) gsub("`", "", x)), 
                                  function(x) nlevels(factor(data[[x]])) - 1))
  }
  if (sum(design.type == 4) > 0) {
    npar.o <- npar.o + sum(sapply(sapply(termlabels[design.type == 4], function(x) gsub("`", "", x)),
                                  function(x) 2^nlevels(factor(data[[x]])) - 2))
    npar.n <- npar.n + sum(sapply(sapply(termlabels[design.type == 4], function(x) gsub("`", "", x)),
                                  function(x) nlevels(factor(data[[x]])) - 1))
  }

  coefs.int <- coefs.int.reparam <- coefs.boundary <- nn.coefs[1]
  coefs.model <- nn.coefs[-1L]
  coefs.model.reparam <- rep(0, npar.n - 1)
  coef.count.o <- coef.count.n <- coef.cont.count <- 0L
    
  coef.names <- "(Intercept)"
    
  for (i in seq_len(length(design.ref))) {
    varname <- gsub("`", "", termlabels[i])
    varref <- allref[[termlabels[i]]][[ref.vector[i]]]
    if (design.type[i] == 1) {
      thiscoef <- coefs.model[(coef.count.o + 1L)]
      cont.min <- min(data[[varname]])
      cont.max <- max(data[[varname]])
      if (type == "cem" || mono[termlabels[i]]) {
        if (mono[termlabels[i]]) coefs.boundary <- append(coefs.boundary, thiscoef)
        if (varref == 1) {
          coefs.int.reparam <- coefs.int.reparam - thiscoef * cont.min
          coefs.model.reparam[(coef.count.n + 1L)] <- thiscoef
        } else {
          coefs.int.reparam <- coefs.int.reparam + thiscoef * cont.max
          coefs.model.reparam[(coef.count.n + 1L)] <- -thiscoef
        }
      } else {
        thiscoef.rev <- coefs.model[(npar.o + coef.cont.count)]
        coefs.int.reparam <- coefs.int.reparam - thiscoef * cont.min + thiscoef.rev * cont.max
        coefs.model.reparam[(coef.count.n + 1L)] <- thiscoef - thiscoef.rev
      }
      coef.names <- append(coef.names, varname)
      coef.count.o <- coef.count.o + 1L
      coef.count.n <- coef.count.n + 1L
      coef.cont.count <- coef.cont.count + as.numeric(!mono[termlabels[i]])
    } else if (design.type[i] == 2) {
      lev <- levels(factor(data[[varname]]))
      nlev <- nlevels(factor(data[[varname]]))
      ref.orig <- varref
      ref.new <- lev[1L]
      if (type == "cem") {
        if (ref.orig != ref.new) {
          thiscoef <- coefs.model[(coef.count.o + 1L):(coef.count.o + nlev - 1L)]
          thiscoef.new <- append(thiscoef, 0, after = which(lev == ref.orig) - 1L)[-1L] - thiscoef[1L]
          coefs.int.reparam <- coefs.int.reparam + thiscoef[1L]
          coefs.model.reparam[(coef.count.n + 1L):(coef.count.n + nlev - 1L)] <- thiscoef.new
        }
      } else {
        thiscoef <- coefs.model[(coef.count.o + 1L):(coef.count.o + nlev)]
        thiscoef.new <- thiscoef[-1L] - thiscoef[1L]
        coefs.int.reparam <- coefs.int.reparam + thiscoef[1L]
        coefs.boundary[1L] <- coefs.boundary[1L] + min(thiscoef)
        coefs.model.reparam[(coef.count.n + 1L):(coef.count.n + nlev - 1L)] <- thiscoef.new
      }
      coef.names <- append(coef.names, paste0(varname, lev[-1L]))
      coef.count.o <- coef.count.o + nlev - as.numeric(type == "cem")
      coef.count.n <- coef.count.n + nlev - 1L
    } else if (design.type[i] == 3) {
      lev <- levels(factor(data[[varname]]))
      nlev <- nlevels(factor(data[[varname]]))
      ref.orig <- varref
      thiscoef <- coefs.model[(coef.count.o + 1L):(coef.count.o + nlev - 1L)]
      thiscoef.sum <- append(cumsum(thiscoef), 0, after = 0)
      thiscoef.ord <- thiscoef.sum[match(lev, ref.orig)]
      thiscoef.new <- thiscoef.ord[-1L] - thiscoef.ord[1L]
      coefs.int.reparam <- coefs.int.reparam + thiscoef.ord[1L]
      coefs.model.reparam[(coef.count.n + 1L):(coef.count.n + nlev - 1L)] <- thiscoef.new
      coef.names <- append(coef.names, paste0(varname, lev[-1L]))
      coef.count.o <- coef.count.o + nlev - 1L
      coef.count.n <- coef.count.n + nlev - 1L
    } else if (design.type[i] == 4) {
      lev <- levels(factor(data[[varname]]))
      nlev <- nlevels(factor(data[[varname]]))
      ref.orig <- varref
      allcombin <- lapply(seq_len(nlev - 1), function(x) combinat::combn(ref.orig, x, simplify = FALSE))
      allcombin2 <- unlist(allcombin, recursive = FALSE)
      thiscoef <- coefs.model[(coef.count.o + 1L):(coef.count.o + length(allcombin2))]
      reord <- rev(cumsum(sapply(allcombin, length))) - 1
      thiscoef.reord <- thiscoef
      for (j in seq_along(reord)) 
        thiscoef.reord <- append(thiscoef.reord[-1L], thiscoef.reord[1L], after = reord[j])
      thiscoef.all <- rep(0, nlev)
      for (k in seq_along(lev))
        thiscoef.all[k] <- sum(thiscoef.reord[sapply(allcombin2, function(x) lev[k] %in% x)])
      thiscoef.new <- thiscoef.all[-1L] - thiscoef.all[1L]
      coefs.int.reparam <- coefs.int.reparam + thiscoef.all[1L]
      coefs.boundary[1L] <- coefs.boundary[1L] + min(thiscoef.all)
      coefs.model.reparam[(coef.count.n + 1L):(coef.count.n + nlev - 1L)] <- thiscoef.new
      coef.names <- append(coef.names, paste0(varname, lev[-1L]))
      coef.count.o <- coef.count.o + length(allcombin2)
      coef.count.n <- coef.count.n + nlev - 1L
    }
  }
  coefs <- c(coefs.int.reparam, coefs.model.reparam)
  names(coefs) <- coef.names
    
  design <- model.matrix(terms, data) 
  list(coefs = coefs, design = design, coefs.boundary = coefs.boundary)
}