R/conTest_ceq_conGLM_F_Wald_score.R

In LeonardV/restriktor: Restricted Statistical Estimation and Inference for Linear Models

conTest_ceq.conGLM <- function(object, test = "F", boot = "no", 
                               R = 9999, p.distr = rnorm, 
                               parallel = "no", ncpus = 1L, cl = NULL, 
                               seed = 1234, verbose = FALSE, ...) {
  
  if (!inherits(object, "conGLM")) {
    stop("object must be of class conGLM.")
  }
  
  test <- tolower(test)
  stopifnot(test %in% c("f","wald","score"))
  
  CON  <- object$CON
  CON$Amat <- Amat <- object$constraints
  CON$bvec <- object$rhs
  CON$meq  <- meq <- object$neq
  
  if (#length(CON$ceq.linear.idx)     > 0  && # some linear eq. constraints
    #length(CON$ceq.nonlinear.idx) == 0L && # no nonlinear eq. constraints
    #length(CON$cin.linear.idx)    == 0L && # no inequality constraints
    #length(CON$cin.nonlinear.idx) == 0L
    nrow(Amat) == meq) {
    
    if (test == "default") {
      test <- "f"
    }
    
    # here we perform the usual Wald/F test...
    if (test == "wald") {
      #theta.r <- object$b.unrestr
      Wald.out <- con_test_Wald(Sigma   = object$Sigma,
                                JAC     = Amat,         
                                theta.r = Amat %*% object$b.unrestr - object$rhs) 
      
      OUT <- append(CON, Wald.out)
      OUT$b.restr <- object$b.restr
      OUT$b.unrestr <- object$b.unrestr
    } else if (test == "f") {
      Wald.out <- con_test_Wald(Sigma   = object$Sigma,
                                JAC     = Amat,
                                theta.r = Amat%*%object$b.unrestr - object$rhs) 
      # convert Wald to F
      OUT <- append(CON, Wald.out)
      OUT$test <- "F"
      OUT$Ts <- Wald.out$Ts / Wald.out$df
      OUT$df.residual <- object$df.residual
      OUT$pvalue <- 1 - pf(OUT$Ts, OUT$df, OUT$df.residual)
      OUT$b.restr <- object$b.restr
      OUT$b.unrestr <- object$b.unrestr
    } else if (test == "score") {
      OUT <- CON
      OUT$test <- "Score"
      # model matrix
      X <- model.matrix(object)[,,drop = FALSE]
      #n <- dim(X)[1]
      #p <- dim(X)[2]
      
      # residuals under the null-hypothesis
      res0 <- residuals(object, "working")
      wres0 <- as.vector(res0) * weights(object, "working")
      # degrees-of-freedom under the null-hypothesis
      #df0 <- (n - (p - qr(Amat[0:meq,,drop = FALSE])$rank))
      # dispersion parameter under the null-hypothesis
      dispersion0 <- object$dispersion
      
      # information matrix under the null-hypothesis
      I0 <- object$information
      # score vector
      G0 <- colSums(as.vector(wres0) * X) / dispersion0
      # score test-statistic
      OUT$Ts <- G0 %*% solve(I0, G0)
      # df
      OUT$df <- nrow(Amat)
      # p-value based on chisq
      OUT$pvalue    <- 1 - pchisq(OUT$Ts, df = OUT$df)
      OUT$b.restr   <- object$b.restr
      OUT$b.unrestr <- object$b.unrestr
    } else {
      stop("restriktor ERROR: test ", sQuote(test), " not (yet) implemented.")
    }
  } else if (#length(CON$ceq.nonlinear.idx) == 0L &&
    #length(CON$cin.linear.idx)     > 0L && # some inequalities restr.
    #length(CON$cin.nonlinear.idx) == 0L
    nrow(Amat != meq)) {
    stop("test not applicable with inequality constraints.")
  } else if (length(CON$ceq.nonlinear.idx) > 0L || length(CON$cin.nonlinear.idx) > 0L) {
    stop("ERROR: can not handle (yet) nonlinear (in)equality constraints")
  }
  
  OUT$boot <- boot
  
  if (boot == "parametric") {
    
    if (!is.function(p.distr)) {
      p.distr <- get(p.distr, mode = "function")
    }
    arguments <- list(...)
    pnames <- names(formals(p.distr))
    pm <- pmatch(names(arguments), pnames, nomatch = 0L)
    pm <- names(arguments)[pm > 0L]
    formals(p.distr)[pm] <- unlist(arguments[pm])
    
    OUT$pvalue <- con_pvalue_boot_parametric(object,
                                             Ts.org   = OUT$Ts,
                                             type     = "A",
                                             test     = test,
                                             R        = R,
                                             p.distr  = p.distr,
                                             parallel = parallel,
                                             ncpus    = ncpus,
                                             cl       = cl,
                                             seed     = seed,
                                             verbose  = verbose)
  } else if (boot == "model.based") {
    stop("Restriktor ERROR: the model-based bootstrap is not implemented.")
  } 
  
  OUT$model.org <- object$model.org
  
  class(OUT) <- "conTest"
  
  OUT
}



con_test_Wald <- function(Sigma, JAC, theta.r) {
  
  # remove redundant rows from JAC *and* theta_r
  npar <- ncol(JAC)
  
  JAC.aug <- cbind(JAC, theta.r)
  Q <- qr(t(JAC.aug))
  JAC.full <- t(qr.X(Q, ncol = Q$rank))
  JAC <- JAC.full[,seq_len(npar),drop = FALSE]
  theta.r <- as.numeric(JAC.full[,(npar + 1L)])
  
  # restricted vcov
  info.r  <- JAC %*% Sigma %*% t(JAC)
  
  # Wald test statistic
  Wald <- as.numeric(t(theta.r) %*% solve(info.r) %*% theta.r)
  
  # df
  Wald.df <- nrow(JAC)
  
  # p-value based on chisq
  Wald.pvalue <- 1 - pchisq(Wald, df = Wald.df)
  
  OUT <- list(test   = "Wald",
              Ts     = Wald,
              df     = Wald.df,
              pvalue = Wald.pvalue)
  
  OUT
}