R/step3_1t.R

In LTCDM: Latent Transition Cognitive Diagnosis Model with Covariates

#==================The modified objective function=====================
obj.logitReg=function(par, w, Z){ #Z: covariate(N*(C+1))； par: initial value of regression coefficient
  #Z already has intercept column
  P=exp(Z%*%par)/(1+exp(Z%*%par))
  ll=-sum( log( (1-P)*w[,1] + P*w[,2] ) )
  return(ll)
}



#' @title Result function for latent regression
#' @description
#' Result function for latent regression.Details can be found in Iaconangelo (2017).
#' @references Iaconangelo, C. (2017). Uses of classification error probabilities in the three-step approach to estimating cognitive diagnosis models (Doctoral dissertation). https://rucore.libraries.rutgers.edu/rutgers-lib/55495/PDF/1/play/
#'
#' @param eap expected a posteriori.
#' @param cep classification error probabilities.
#' @param covar covariances.
#' @param K number of attributes.
#' @param method the sample-level \code{"SL"} or the posterior distribution level \code{"PDL"} correction weights.
#'
#' @return A data.frame containing the following columns:
#' \describe{
#'   \item{Attribute}{The attribute.}
#'   \item{Estimate}{Estimated regression coefficient.}
#'   \item{odds.ratio}{Odds ratio.}
#'   \item{d}{Effect size.}
#'   \item{SE}{Standard error.}
#'   \item{CI.95}{95% confidence interval.}
#'   \item{z.value}{Wald test z-statistics.}
#'   \item{p.2tailed}{Two-tailed p-value.}
#'   \item{p.1tailed}{One-tailed p-value.}
#' }
#'
#' @export
#' @examples
#' # ---- Runable example (using simulated CEP) ----
#' if (requireNamespace("GDINA", quietly = TRUE)) {
#'   library(GDINA)
#'   dat <- sim10GDINA$simdat
#'   Q <- matrix(c(1,0,0,
#'                 0,1,0,
#'                 0,0,1,
#'                 1,0,1,
#'                 0,1,1,
#'                 1,1,0,
#'                 0,0,1,
#'                 1,0,0,
#'                 1,1,1,
#'                 1,0,1), byrow = TRUE, ncol = 3)
#'   fit.object <- GDINA(dat, Q = Q, model = "GDINA", att.dist = "independent", verbose = FALSE)
#'   cep_test <- list(SL.k = list(), PDL.k = list())
#'   for (i in 1:3) {
#'     cep_test$SL.k[[i]] <- matrix(runif(4, min = 0, max = 1), nrow = 2)
#'     cep_test$PDL.k[[i]] <- lapply(1:1000, function(j) matrix(runif(4, min = 0, max = 1), nrow = 2))
#'   }
#'   eap <- personparm(fit.object, what = "EAP")
#'   covar_test <- matrix(runif(3000, min = -2, max = 2), nrow = 1000)
#'   K <- 3
#'   result_SL <- step3_1t(eap = eap, cep = cep_test, covar = covar_test, K = K, method = "SL")
#'   result_PDL <- step3_1t(eap = eap, cep = cep_test, covar = covar_test, K = K, method = "PDL")
#' }
#'
#' # ---- Not recommended to run (depends on CEP() output) ----
#' \dontrun{
#' fit.object <- GDINA(dat = dat, Q = Q, model = "GDINA", att.dist = "independent", verbose = FALSE)
#' cep <- CEP_1t(fit.object)
#' eap <- personparm(fit.object, what = "EAP")
#' out_PDL <- step3_1t(eap, Q, cep, covar, K, "PDL")
#' out_SL <- step3_1t(eap, Q, cep, covar, K, "SL")
#' }
#'



step3_1t <- function(eap,cep,covar,K,method){

  C = ncol(covar) #number of covariates
  N = nrow(covar)
  B = matrix(1,C,K) #initial values of regression coefficients
  res = list()
  res.k = list()

  for (k in 1:K) {

    ind=c(eap[,k]+1)
    w.SL=t(cep$SL.k[[k]][,ind])

    w.PDL=matrix(NA,N,2)
    for (i in 1:N){
      w.PDL[i,]=cep$PDL.k[[k]][[i]][,ind[i]]
    }

    #PDL weight should be better
    #================= Optimize the modified objective function =================
    out.SL=optim(par=as.matrix(B[,k]),fn=obj.logitReg, w=w.SL, Z=covar,method="L-BFGS-B",
                 hessian=T, lower=-5, upper=5)

    out.PDL=optim(par=as.matrix(B[,k]),fn=obj.logitReg, w=w.PDL, Z=covar, method="L-BFGS-B",
                  hessian=T,lower=-5, upper=5)


    if(method == "SL"){
      out.cor = out.SL
    }else if(method == "PDL"){
      out.cor = out.PDL
    }

    #compute SEs of coefficients
    mathessian <- out.cor$hessian
    inversemathessian <- solve(mathessian)
    SE <- sqrt(diag(inversemathessian)) #check this

    # Wald test
    z.value <- out.cor$par / SE
    p.1tailed <- pnorm(-abs(z.value))
    p.2tailed <- 2 * pnorm(-abs(z.value))

    # Confidence Intervals
    CI_l <- round(out.cor$par - 1.96 * SE, 4)
    CI_r <- round(out.cor$par + 1.96 * SE, 4)
    CI.95 <- paste0("(", CI_l, ",", CI_r, ")")

    # results
    # Compile result
    res[[k]] <- data.frame(
      Attribute = rep(k, length(out.cor$par)),
      Estimate = round(out.cor$par, 4),
      odds.ratio = round(exp(out.cor$par), 4),
      d = round(out.cor$par * sqrt(3)/pi, 4),
      SE = round(SE, 4),
      CI.95 = CI.95,
      z.value = round(z.value, 4),
      p.2tailed = round(p.2tailed, 4),
      p.1tailed = round(p.1tailed, 4)
    )
  }
    result <- do.call(rbind, res)

    return(result)

}