R/help_fit.R

In GPCMlasso: Regularized Explanatory Generalized Partial Credit Models

loglikscoreDIFlasso2 <- function(alpha, Y, X, Z, Q, q, n, I,
                                 px, GHweights, GHnodes,
                                 acoefs, lambda, lambda2, cvalue, cores,
                                 weight, n_sigma, scale_fac,
                                 scale_cols = numeric(0)) {
  
  ## scale_cols wird hier ignoriert, damit generische Aufrufe nicht scheitern.
  l <- loglikscoreDIFlasso(alpha, Y, X, Z, Q, q, n, I,
                           px, GHweights, GHnodes,
                           acoefs, lambda, lambda2, cvalue, cores,
                           weight, n_sigma, scale_fac)
  
  ret <- l$objective
  attr(ret, "gradient") <- l$gradient
  ret
}



loglikscorePCMlasso2 <- function(alpha, Y, X, Z, Q, q, n, I,
                                 px, GHweights, GHnodes,
                                 acoefs, lambda, lambda2, cvalue, cores,
                                 weight, n_sigma, scale_fac,
                                 scale_cols = numeric(0)) {
  
  l <- loglikscorePCMlasso(alpha, Y, X, Z, Q, q, n, I,
                           px, GHweights, GHnodes,
                           acoefs, lambda, lambda2, cvalue, cores,
                           weight, n_sigma, scale_fac,
                           scale_cols = scale_cols)
  
  ret <- l$objective
  attr(ret, "gradient") <- l$gradient
  ret
}



## ------------------------------------------------------------
## Helper: scale_cols passend zu aktueller Z-Matrix machen
##
## combined design in C++ ist immer cbind(design, Z)
##
## scale_cols bezieht sich aber auf cbind(design, designX)
## und muss deshalb angepasst werden, wenn statt designX z.B.
## design_null verwendet wird.
## ------------------------------------------------------------
make_scale_cols_current <- function(scale_cols, design, designX, Z_current) {
  
  if (length(scale_cols) == 0) {
    return(numeric(0))
  }
  
  n_design <- ncol(design)
  n_Z_full <- ncol(designX)
  n_Z_current <- ncol(Z_current)
  
  if (is.null(n_Z_current)) {
    n_Z_current <- 0
  }
  
  if (length(scale_cols) != n_design + n_Z_full) {
    stop("scale_cols must have length ncol(design) + ncol(designX).")
  }
  
  scale_design <- scale_cols[seq_len(n_design)]
  
  if (n_Z_current == 0) {
    return(scale_design)
  }
  
  ## Assumption: Z_current contains the first n_Z_current columns of designX
  scale_Z <- scale_cols[n_design + seq_len(n_Z_full)]
  
  return(c(scale_design, scale_Z[seq_len(n_Z_current)]))
}



## ------------------------------------------------------------
## Helper: Check, ob eine Funktion scale_cols als Argument hat
## ------------------------------------------------------------
fun_has_scale_cols <- function(fun) {
  "scale_cols" %in% names(formals(fun))
}



help_fit <- function(model, Y, l.lambda, start, loglik_fun, score_fun, log_score_fun, adaptive,
                     Q, q, I, n, m, response, design, designX, px,
                     GHweights, GHnodes, acoefs, lambda2, cvalue, n_sigma,
                     l.bound, trace, log.lambda, weight.penalties, scale_fac = scale_fac,
                     ada.lambda, lambda.min, ada.power, cores,
                     null_thresh, DSF, gradtol, iterlim, steptol,
                     main.effects, penalize.main.effects, ctrl.gpcm,
                     scale_cols = numeric(0)) {
  
  ## get initial weight parameters
  weight <- rep(1, ncol(acoefs))
  
  
  ## initialize starting values
  if (is.null(start)) {
    
    if (trace) {
      cat("Find start values ...\n")
    }
    
    alpha.start <- c(rep(0.1, px))
    
    
    if (!(model %in% c("RSM", "GRSM"))) {
      
      if (model == "GPCM") {
        m.ltm <- gpcm(Y, constraint = "gpcm", control = ctrl.gpcm)
        coef.ltm <- coef(m.ltm)
        
        if (is.matrix(coef.ltm)) {
          sigma.start <- coef.ltm[, q[1] + 1]
          delta.start <- c(t(coef.ltm[, -(q[1] + 1)]))
        } else {
          sigma.start <- delta.start <- c()
          for (u in 1:I) {
            sigma.start[u] <- coef.ltm[[u]][q[u] + 1]
            delta.start <- c(delta.start, coef.ltm[[u]][-(q[u] + 1)])
          }
        }
      }
      
      
      if (model == "PCM") {
        m.ltm <- gpcm(Y, constraint = "1PL", control = ctrl.gpcm)
        coef.ltm <- coef(m.ltm)
        
        if (is.matrix(coef.ltm)) {
          sigma.start <- coef.ltm[1, q[1] + 1]
          delta.start <- c(t(coef.ltm[, -(q[1] + 1)]))
        } else {
          delta.start <- c()
          sigma.start <- coef.ltm[[1]][q[1] + 1]
          for (u in 1:I) {
            delta.start <- c(delta.start, coef.ltm[[u]][-(q[u] + 1)])
          }
        }
      }
      
      
      if (model == "2PL") {
        m.ltm <- gpcm(Y, constraint = "gpcm", control = ctrl.gpcm)
        coef.ltm <- coef(m.ltm)
        sigma.start <- coef.ltm[, q[1] + 1]
        delta.start <- c(t(coef.ltm[, -(q[1] + 1)]))
      }
      
      
      if (model == "RM") {
        m.ltm <- gpcm(Y, constraint = "1PL", control = ctrl.gpcm)
        coef.ltm <- coef(m.ltm)
        sigma.start <- coef.ltm[1, q[1] + 1]
        delta.start <- c(t(coef.ltm[, -(q[1] + 1)]))
      }
      
    } else {
      
      if (model == "RSM") {
        m.mirt <- mirt(Y, 1, itemtype = "rsm", verbose = FALSE)
        coefmethod <- selectMethod("coef", class(m.mirt))
        coef.mirt <- coefmethod(m.mirt, simplify = TRUE)
        sigma.start <- sqrt(coef.mirt$cov)
        alpha.mirt <- (coef.mirt$items)[1, 2:(q[1] + 1)]
        delta.mirt <- (coef.mirt$items)[, (q[1] + 2)]
        delta.mirt <- -delta.mirt + alpha.mirt[1]
        alpha.mirt <- alpha.mirt - alpha.mirt[1]
        delta.start <- c(delta.mirt, alpha.mirt[-1])
      }
      
      
      if (model == "GRSM") {
        m.mirt <- mirt(Y, 1, itemtype = "rsm", verbose = FALSE)
        coefmethod <- selectMethod("coef", class(m.mirt))
        coef.mirt <- coefmethod(m.mirt, simplify = TRUE)
        sigma.start <- rep(sqrt(coef.mirt$cov), n_sigma)
        alpha.mirt <- (coef.mirt$items)[1, 2:(q[1] + 1)]
        delta.mirt <- (coef.mirt$items)[, (q[1] + 2)]
        delta.mirt <- -delta.mirt + alpha.mirt[1]
        alpha.mirt <- alpha.mirt - alpha.mirt[1]
        delta.start <- c(delta.mirt, alpha.mirt[-1])
      }
    }
    
    
    alpha.start <- c(delta.start, rep(0, ncol(designX)), abs(sigma.start))
    
    alpha.null <- alpha.start[rowSums(abs(acoefs)) == 0]
    p_null <- length(alpha.null)
    
    design_null <- matrix(0, 0, 0)
    
    if (main.effects & ncol(designX) > 0 & (!penalize.main.effects)) {
      design_null <- designX[, 1:m, drop = FALSE]
    }
    
    scale_cols_null <- make_scale_cols_current(
      scale_cols = scale_cols,
      design = design,
      designX = designX,
      Z_current = design_null
    )
    
    acoefs_null <- matrix(0, nrow = p_null, ncol = 1)
    bound_null <- l.bound[rowSums(abs(acoefs)) == 0]
    
    loglik_NA <- TRUE
    
    
    while (loglik_NA) {
      
      loglik_NA <- is.nan(
        log_score_fun(alpha.null,
                      Q = Q, q = q, I = I, n = n,
                      Y = response,
                      X = design,
                      Z = design_null,
                      px = p_null,
                      GHweights = GHweights,
                      GHnodes = GHnodes,
                      acoefs = acoefs_null,
                      lambda = 0,
                      scale_fac = scale_fac,
                      lambda2 = lambda2,
                      cvalue = cvalue,
                      cores = cores,
                      weight = 1,
                      n_sigma = n_sigma,
                      scale_cols = scale_cols_null)
      )
      
      if (loglik_NA) {
        alpha.null <- alpha.null * 0.9
      }
    }
    
    
    if (fun_has_scale_cols(loglik_fun)) {
      
      m.opt <- try(
        nlminb(start = alpha.null,
               objective = loglik_fun,
               gradient = score_fun,
               Q = Q,
               q = q,
               I = I,
               n = n,
               Y = response,
               X = design,
               Z = design_null,
               px = p_null,
               GHweights = GHweights,
               GHnodes = GHnodes,
               acoefs = acoefs_null,
               lambda = 0,
               scale_fac = scale_fac,
               lambda2 = lambda2,
               cvalue = cvalue,
               cores = cores,
               weight = 1,
               n_sigma = n_sigma,
               scale_cols = scale_cols_null,
               control = list(eval.max = 500, iter.max = 500, step.min = 0.01),
               lower = bound_null)
      )
      
    } else {
      
      m.opt <- try(
        nlminb(start = alpha.null,
               objective = loglik_fun,
               gradient = score_fun,
               Q = Q,
               q = q,
               I = I,
               n = n,
               Y = response,
               X = design,
               Z = design_null,
               px = p_null,
               GHweights = GHweights,
               GHnodes = GHnodes,
               acoefs = acoefs_null,
               lambda = 0,
               scale_fac = scale_fac,
               lambda2 = lambda2,
               cvalue = cvalue,
               cores = cores,
               weight = 1,
               n_sigma = n_sigma,
               control = list(eval.max = 500, iter.max = 500, step.min = 0.01),
               lower = bound_null)
      )
    }
    
    
    if (inherits(m.opt, "try-error")) {
      stop("Initial null model estimation failed.")
    }
    
    
    alpha.start[rowSums(abs(acoefs)) == 0] <- m.opt$par
    alpha.start[rowSums(abs(acoefs)) != 0] <- 1e-8
    
  } else {
    
    alpha.start <- start
  }
  
  
  ## get new weights if necessary
  if (adaptive) {
    
    if (trace) {
      cat("Get adaptive weights ...", "\n")
    }
    
    scale_cols_full <- make_scale_cols_current(
      scale_cols = scale_cols,
      design = design,
      designX = designX,
      Z_current = designX
    )
    
    
    if (fun_has_scale_cols(loglik_fun)) {
      
      m.opt <- try(
        nlminb(start = alpha.start,
               objective = loglik_fun,
               gradient = score_fun,
               Q = Q,
               q = q,
               I = I,
               n = n,
               Y = response,
               X = design,
               Z = designX,
               px = px,
               GHweights = GHweights,
               GHnodes = GHnodes,
               acoefs = acoefs,
               lambda = 0,
               scale_fac = scale_fac,
               lambda2 = ada.lambda,
               cvalue = cvalue,
               cores = cores,
               weight = weight,
               n_sigma = n_sigma,
               scale_cols = scale_cols_full,
               control = list(eval.max = 500, iter.max = 500, step.min = 0.01),
               lower = l.bound)
      )
      
    } else {
      
      m.opt <- try(
        nlminb(start = alpha.start,
               objective = loglik_fun,
               gradient = score_fun,
               Q = Q,
               q = q,
               I = I,
               n = n,
               Y = response,
               X = design,
               Z = designX,
               px = px,
               GHweights = GHweights,
               GHnodes = GHnodes,
               acoefs = acoefs,
               lambda = 0,
               scale_fac = scale_fac,
               lambda2 = ada.lambda,
               cvalue = cvalue,
               cores = cores,
               weight = weight,
               n_sigma = n_sigma,
               control = list(eval.max = 500, iter.max = 500, step.min = 0.01),
               lower = l.bound)
      )
    }
    
    
    if (inherits(m.opt, "try-error")) {
      stop("Adaptive weights can not be calculated! Increase ada.lambda or set adaptive = FALSE!")
    }
    
    
    weight <- try(m.opt$par)
    weight <- abs(t(acoefs) %*% weight)^ada.power
    
    if (any(weight == 0)) {
      weight[which(weight == 0)] <- 1e-4
    }
    
    weight <- as.vector(1 / weight)
  }
  
  
  ## find maximal lambda value and make grid
  if (!is.na(l.lambda)) {
    
    if (trace) {
      cat("Find maximal tuning parameter ...", "\n")
    }
    
    scale_cols_full <- make_scale_cols_current(
      scale_cols = scale_cols,
      design = design,
      designX = designX,
      Z_current = designX
    )
    
    
    if (fun_has_scale_cols(score_fun)) {
      
      score <- score_fun(alpha.start,
                         response,
                         design,
                         designX,
                         Q,
                         q,
                         n,
                         I,
                         px,
                         GHweights,
                         GHnodes,
                         acoefs,
                         0,
                         lambda2,
                         cvalue,
                         cores,
                         weight,
                         n_sigma,
                         scale_fac,
                         scale_cols = scale_cols_full)
      
    } else {
      
      score <- score_fun(alpha.start,
                         response,
                         design,
                         designX,
                         Q,
                         q,
                         n,
                         I,
                         px,
                         GHweights,
                         GHnodes,
                         acoefs,
                         0,
                         lambda2,
                         cvalue,
                         cores,
                         weight,
                         n_sigma,
                         scale_fac)
    }
    
    
    a <- abs(score / acoefs %*% weight)
    a[a == Inf] <- 0
    lambda.max <- max(a[rowSums(abs(acoefs)) != 0]) * 1.1
    
    
    if (DSF) {
      
      ## Falls DSF = TRUE verwendet wird, muss find.lambda separat
      ## an scale_cols angepasst werden.
      ## Fuer den aktuellen GPCM-Test ohne DSF ist dieser Zweig irrelevant.
      lambda.max <- find.lambda(lambda.max,
                                l.lambda,
                                alpha.start,
                                log_score_fun,
                                Q,
                                q,
                                I,
                                n,
                                response,
                                design,
                                designX,
                                px, 
                                GHweights,
                                GHnodes,
                                acoefs,
                                scale_fac,
                                lambda2,
                                cvalue, 
                                cores,
                                weight,
                                n_sigma,
                                null_thresh,
                                gradtol,
                                iterlim, 
                                steptol)
    }
    
    
    if (log.lambda) {
      
      correct.factor <- 0.0
      
      lambda <- exp(seq(log(lambda.max + correct.factor * lambda.max), 
                        log(lambda.min + correct.factor * lambda.max),
                        length = l.lambda)) -
        correct.factor * lambda.max
      
      lambda[l.lambda] <- lambda.min
      
    } else {
      
      lambda <- seq(lambda.max, lambda.min, length = l.lambda)
    }
    
  } else {
    
    lambda <- NA
  }
  
  
  return(list(lambda = lambda,
              weight = weight,
              alpha.start = alpha.start))
}