psindex: Parameter Stability Index

# This code is contributed by Leonard Vanbrabant <L.G.F.Vanbrabant@gmail.com>
InformativeTesting <- function(model = NULL, data, constraints = NULL, 
                               R = 1000L, type = "bollen.stine",
                               return.LRT = TRUE, 
                               double.bootstrap = "standard",
                               double.bootstrap.R = 249,
                               double.bootstrap.alpha = 0.05,
                               parallel = c("no", "multicore", "snow"), 
                               ncpus = 1L, cl = NULL, verbose = FALSE, ...){
  
  fit.B1 <- sem(model, ...,
                data = data, 
                se   = "none", 
                test = "standard") 
  
  fit.B0 <- fit.A1 <- sem(model, ...,
                          data        = data, 
                          se          = "none", 
                          test        = "standard", 
                          constraints = constraints) 
  
  con.idx  <- (max(fit.B1@ParTable$id) + 1L):max(fit.A1@ParTable$id)
  
  user.equal  <- fit.A1@ParTable
  user.equal$op[con.idx] <- "=="
  
  fit.A0 <- sem(user.equal, ...,
                data = data, 
                se   = "none", 
                test = "standard")
  
  lrt.bootA <- bootstrapLRT(fit.A0, fit.A1, 
                            R                      = R, 
                            type                   = type, 
                            verbose                = verbose,
                            return.LRT             = return.LRT, 
                            double.bootstrap       = double.bootstrap,
                            double.bootstrap.R     = double.bootstrap.R, 
                            double.bootstrap.alpha = double.bootstrap.alpha,
                            parallel               = parallel, 
                            ncpus                  = ncpus, 
                            cl                     = cl)
  
  lrt.bootB <- bootstrapLRT(fit.B0, fit.B1, 
                            R                      = R, 
                            type                   = type, 
                            verbose                = verbose,
                            return.LRT             = return.LRT, 
                            double.bootstrap       = double.bootstrap,
                            double.bootstrap.R     = double.bootstrap.R, 
                            double.bootstrap.alpha = double.bootstrap.alpha,
                            parallel               = parallel, 
                            ncpus                  = ncpus, 
                            cl                     = cl)
  
  output <- list(fit.A0 = fit.A0, fit.A1 = fit.A1, fit.B1 = fit.B1,
                 lrt.bootA = lrt.bootA, lrt.bootB = lrt.bootB,
                 double.bootstrap = double.bootstrap,
                 double.bootstrap.alpha = double.bootstrap.alpha,
                 return.LRT = return.LRT, type = type)
  
  class(output) <- "InformativeTesting"
  
  return(output)
}


print.InformativeTesting <- function(x, digits = max(3, getOption("digits") - 3), ...) {
  object <- x
  cat("\nInformativeTesting: Order/Inequality Constrained Hypothesis Testing:\n\n")
  cat("  Variable names in model         :", unlist(object$fit.A1@Data@ov.names[1]), "\n")  
  cat("  Number of variables             :", object$fit.A1@Model@nvar[1], "\n")  
  cat("  Number of groups                :", object$fit.A1@Data@ngroups, "\n")  
  cat("  Used sample size per group      :", unlist(object$fit.A1@Data@nobs), "\n")
  cat("  Used sample size                :", sum(unlist(object$fit.A1@Data@nobs)), "\n")
  cat("  Total sample size               :", sum(unlist(object$fit.A1@Data@norig)), "\n\n")
  cat("  Estimator                       :", object$fit.A1@Options$estimator, "\n")
  cat("  Missing data                    :", object$fit.A1@Options$missing, "\n")
  cat("  Bootstrap method                :", object$type, "\n")
  cat("  Double bootstrap method         :", object$double.bootstrap, "\n")
  
  dbtype <- object$double.bootstrap
  # original LRT for hypothesis test Type A
  TsA <- attr(object$lrt.bootA, "LRT.original")
  # original LRT for hypothesis test Type B
  TsB <- attr(object$lrt.bootB, "LRT.original")
  # unadjusted pvalues for Ts
  pvalueA <- object$lrt.bootA[1]
  pvalueB <- object$lrt.bootB[1]
  alpha <- object$double.bootstrap.alpha
  
  ###
  if (dbtype == "no") {
    cat("\n\n  Type A test: H0: all restriktions active (=)", "\n", 
        "          vs. H1: at least one restriktion strictly true (>)", "\n")
    cat("         Test statistic: ", format(round(TsA, digits), nsmall = digits), ", unadjusted p-value: ", 
        if (pvalueA < 1e-04) {
          "<0.0001"
        } else {
          format(round(pvalueA, digits), nsmall = digits)}, " (alpha = ", alpha, ") ", "\n\n", sep = "")
    cat("  Type B test: H0: all restriktions true", "\n", 
        "          vs. H1: at least one restriktion false", "\n")
    cat("         Test statistic: ", format(round(TsB, digits), nsmall = digits), ", unadjusted p-value: ", 
        if (pvalueB < 1e-04) {
          "<0.0001"
        } else {
          format(round(pvalueB, digits), nsmall = digits)}, " (alpha = ", alpha, ") ", "\n", sep = "")
  } else if (dbtype == "FDB") {
    # adjusted pvalues for Ts
    adj.pvalueA <- attr(object$lrt.bootA, "adj.pvalue")
    adj.pvalueB <- attr(object$lrt.bootB, "adj.pvalue")
    cat("\n\n  Type A test: H0: all restriktions active (=)", "\n", 
        "          vs. H1: at least one restriktion strictly true (>)", "\n")
    cat("         Test statistic: ", format(round(TsA, digits), nsmall = digits), ", adjusted p-value: ", 
        if (adj.pvalueA < 1e-04) {
          "<0.0001"
        } else {
          format(round(adj.pvalueA, digits), nsmall = digits)}, " (alpha = ", 
        alpha, ") ", "\n\n", sep = "")
    cat("  Type B test: H0: all restriktions true", "\n", 
        "          vs. H1: at least one restriktion false", "\n")
    cat("         Test statistic: ", format(round(TsB, digits), 
                                            nsmall = digits), ", adjusted p-value: ", 
        if (adj.pvalueB < 1e-04) {
          "<0.0001"
        } else {
          format(round(adj.pvalueB, digits), nsmall = digits)}, " (alpha = ", 
        alpha, ") ", "\n", sep = "")
  } else if (dbtype == "standard") {
    # adjusted nominal levels
    adj.alphaA  <- attr(object$lrt.bootA, "adj.alpha")
    adj.alphaB  <- attr(object$lrt.bootB, "adj.alpha")
    # adjusted pvalues for Ts
    adj.pvalueA <- attr(object$lrt.bootA, "adj.pvalue")
    adj.pvalueB <- attr(object$lrt.bootB, "adj.pvalue")
    cat("\n\n  Type A test: H0: all restriktions active (=)", "\n", 
        "          vs. H1: at least one restriktion strictly true (>)", "\n")
    cat("         Test statistic: ", format(round(TsA, digits), 
                                            nsmall = digits), ", adjusted p-value: ", 
        if (adj.pvalueA < 1e-04) {
          "<0.0001"
        } else {
          format(round(adj.pvalueA, digits), nsmall = digits)}, " (alpha = ", alpha, ") ", "\n", sep = "")
    cat("                                ", "unadjusted p-value: ", 
        if (pvalueA < 1e-04) {
          "<0.0001"
        } else {
          format(round(pvalueA, digits), nsmall = digits)}, " (alpha = ", 
        format(round(adj.alphaA, digits), nsmall = digits), ") ", "\n\n", sep = "")
    cat("  Type B test: H0: all restriktions true", "\n", 
        "          vs. H1: at least one restriktion false", "\n")
    cat("         Test statistic: ", format(round(TsB, digits), nsmall = digits), ", adjusted p-value: ", 
        if (adj.pvalueB < 1e-04) {
          "<0.0001"
        } else {
          format(round(adj.pvalueB, digits), nsmall = digits)}, " (alpha = ", alpha, ") ", "\n", sep = "")
    cat("                               ", "unadjusted p-value: ", 
        if (pvalueB < 1e-04) {
          "<0.0001"
        } else {
          format(round(pvalueB, digits), nsmall = digits)}, " (alpha = ", 
        format(round(adj.alphaB, digits), nsmall = digits), ") ", "\n\n", sep = "")
  }
  
  if (dbtype == "no") {
    cat("\n  No double bootstrap method is set. The results may be spurious.\n\n")
  }
  
}


plot.InformativeTesting <- function(x, ..., 
                                    type       = c("lr", "ppv"),
                                    main       = "main",
                                    xlab       = "xlabel",
                                    ylab       = "Frequency",
                                    freq       = TRUE,
                                    breaks     = 15,
                                    cex.main   = 1,
                                    cex.lab    = 1,
                                    cex.axis   = 1,
                                    col        = "grey",
                                    border     = par("fg"),
                                    vline      = TRUE, 
                                    vline.col  = c("red", "blue"), 
                                    lty        = c(1,2),
                                    lwd        = 1,
                                    legend     = TRUE,
                                    bty        = "o",
                                    cex.legend = 1,
                                    loc.legend = "topright") {
  object <- x
  return.LRT <- object$return.LRT
  double.bootstrap <- object$double.bootstrap
  double.bootstrap.alpha <- object$double.bootstrap.alpha
  pvalue <- c(object$lrt.bootA[1], object$lrt.bootB[1])
  
  par(mfrow = c(1, 2))
  if (length(type) == 2) {
    par(mfrow = c(2, 2))
  }
  
  if (return.LRT && (type == "lr" || length(type) == 2)) {
    lrt.obs  <- c(attr(object$lrt.bootA, "LRT.original"), 
                  attr(object$lrt.bootB, "LRT.original"))
    lrt.A <- attr(object$lrt.bootA, "LRT")
    lrt.B <- attr(object$lrt.bootB, "LRT")
    if (length(lrt.A) - length(lrt.B) < 0L) {
      lrt <- as.data.frame(cbind(c(lrt.A, rep(as.numeric(NA), length(lrt.B) - 
                                                length(lrt.A))), lrt.B))      
    } else { 
      lrt <- as.data.frame(cbind(lrt.A, c(lrt.B, rep(as.numeric(NA), 
                                                     length(lrt.A) - 
                                                       length(lrt.B)))))
    }
    names(lrt) <- c("lrt.A", " lrt.B")
    
    if (xlab == "xlabel") { 
      xlab.lrt <- c("Bootstrapped LR values")
    }
    if (main == "main") {
      main.lrt <- c("Distr. of LR values - Type A", 
                    "Distr. of LR values - Type B")
    }
    
    for (i in 1:2) {
      plot <- hist(lrt[,i], plot = FALSE, breaks = breaks) 
      plot(plot, ...,
           freq     = freq, 
           main     = main.lrt[i], 
           xlab     = xlab.lrt, 
           ylab     = ylab, 
           cex.axis = cex.axis, 
           cex.main = cex.main, 
           cex.lab  = cex.lab, 
           col      = col, 
           border   = border, 
           axes     = FALSE, 
           xaxt     = "n") 
      
      axis(side = 1)
      axis(side = 2)
      box(lty = 1, col = "black")
      
      if (vline) {
        abline(v   = lrt.obs[i], 
               col = vline.col[1], 
               lty = lty[1], 
               lwd = lwd)
      }  
      if (legend) {
        ppvalue  <- sprintf("%.2f", pvalue[i]) 
        obs.lrt  <- sprintf("%.2f", lrt.obs[i])
        ppval    <- paste0("plug-in p value = ", ppvalue)
        obs.lrt  <- paste0("observed LR = ", obs.lrt)
        legend.obj <- c(obs.lrt, ppval)
        if (!vline) {
          legend(loc.legend, legend.obj, 
                 lty = c(0, 0),   
                 lwd = lwd, 
                 cex = cex.legend, 
                 bty = bty)
        } else {
          legend(loc.legend, legend.obj, 
                 lty = c(lty[1], 0), 
                 col = vline.col[1],  
                 lwd = lwd, 
                 cex = cex.legend, 
                 bty = bty)
        }
      }
    } 
  }
  
  if (double.bootstrap == "standard" && (type == "ppv" || length(type) == 2)) {
    ppvalue.A <- attr(object$lrt.bootA, "plugin.pvalues")
    ppvalue.B <- attr(object$lrt.bootB, "plugin.pvalues")
    adj.a <- c(quantile(ppvalue.A, double.bootstrap.alpha), 
               quantile(ppvalue.B, double.bootstrap.alpha))
    adj.ppv <- c(attr(object$lrt.bootA, "adj.pvalue"), 
                 attr(object$lrt.bootB, "adj.pvalue"))
    if (length(ppvalue.A) - length(ppvalue.B) < 0L) {
      ppv <- as.data.frame(cbind(c(ppvalue.A, rep(NA, length(ppvalue.B) - 
                                                    length(ppvalue.A))), ppvalue.B))
    } else { 
      ppv <- as.data.frame(cbind(ppvalue.A, c(ppvalue.B, rep(NA, length(ppvalue.A) - 
                                                               length(ppvalue.B)))))
    }
    names(ppv) <- c("ppA", "ppB")
    
    if (xlab == "xlabel") {  
      xlab.ppv  <- c("Bootstrapped plug-in p-values")
    }
    if (main == "main") {
      main.ppv  <- c("Distr. of plug-in p-values - Type A", 
                     "Distr. of plug-in p-values - Type B")
    }
    
    for (i in 1:2) {
      plot <- hist(ppv[,i], plot = FALSE, breaks=breaks) 
      plot(plot, ...,
           freq     = freq, 
           main     = main.ppv[i], 
           xlab     = xlab.ppv, 
           ylab     = ylab, 
           cex.axis = cex.axis, 
           cex.main = cex.main, 
           cex.lab  = cex.lab, 
           col      = col, 
           border   = border, 
           axes     = FALSE, 
           xaxt     = "n") 
      
      axis(side = 1, at = seq(0,1,0.1))
      axis(side = 2)
      box(lty = 1, col = "black")
      if (vline) {
        abline(v   = adj.a[i], 
               col = vline.col[1], 
               lty = lty[1], 
               lwd = lwd)
        abline(v   = adj.ppv[i], 
               col = vline.col[2], 
               lty = lty[2], 
               lwd = lwd)
      }
      if (legend) {
        adj.alpha  <- sprintf("%.2f", adj.a[i]) 
        adj.pval   <- sprintf("%.2f", adj.ppv[i])
        adja <- paste0("Adjusted alpha = ", adj.alpha)
        adjp <- paste0("Adjusted p-value = ", adj.pval)
        legend.obj <- c(adja, adjp)
        if (!vline) {
          legend(loc.legend, legend.obj, 
                 lty = 0, 
                 col = vline.col,  
                 lwd = lwd, 
                 cex = cex.legend, 
                 bty = bty)
        } else {
          legend(loc.legend, legend.obj, 
                 lty = lty, 
                 col = vline.col,  
                 lwd = lwd, 
                 cex = cex.legend, 
                 bty = bty)
        }
      }
    }
  }
}