R/prediction.R

In blindedmanuscript/LSBP: Dependent Logit Stick-Breaking Process

#' Predict method for the LSBP
#' 
#' 
#' Predict method for the LSBP estimated using the ECM algorithm.
#' 
#' @param object An object of class \verb{LSBP_EM}.
#' @param type String indicating the type of prediction: \verb{type="mean"},\verb{type="variance"} or \verb{type="cdf"}. See details.
#' @param newdata A new data frame containing the same variables declared in \verb{Formula}. If missing, the dataset provided for estimation is used.
#' @param threshold Only needed if \verb{type="cdf"} is selected. See details.
#' @param ... Further arguments passed to or from other methods.
#' 
#' @details The method \verb{predict.LSBP_ECM} produces predicted values, obtained by evaluating the conditional mean (if \verb{type="mean"}), the conditional variance (if \verb{type="variance"}) or the conditional cumulative distribution function (if \verb{type="cdf"}) at a given \verb{threshold}, after plugging-in the MAP, and using the observations contained in the \verb{newdata} data frame.
#' 
#' @export

predict.LSBP_ECM <- function(object, type="mean", newdata=NULL, threshold=NULL, ...) {
  if (!(type %in% c("mean", "variance","cdf"))) 
    stop("Please provide a valid predictive type")
   if (is.null(newdata)) {
      X1 <- object$data$X1
      X2 <- object$data$X2
   } else {
      f <- object$call
      X1 <- model.matrix(f, data = newdata, rhs = 1)
      if (length(f)[2] == 1) {
         X2 <- X1
      } else {
         X2 <- model.matrix(f, data = newdata, rhs = 2)
      }
   }
   if(type=="mean"){
    if (NCOL(X1) > 1) {
      pred <- pred_mean_multi(X1, X2, object$param$beta_mixing, object$param$beta_kernel)
      
    } else {
        # Univariate case
        pred <- pred_mean(X2, object$param$beta_mixing, object$param$beta_kernel)
    }
   }
  
  if(type=="variance"){
    if (NCOL(X1) > 1) {
      pred <- pred_var_multi(X1, X2, object$param$beta_mixing, object$param$beta_kernel,object$param$tau)
      
    } else {
      # Univariate case
      pred <- pred_var(X2, object$param$beta_mixing, object$param$beta_kernel,object$param$tau)
    }
  }
  
  if(type=="cdf"){
    if (NCOL(X1) > 1) {
      pred <- pred_cdf_multi(X1, X2, object$param$beta_mixing, object$param$beta_kernel,object$param$tau,threshold)
      
    } else {
      # Univariate case
      pred <- pred_cdf(X2, object$param$beta_mixing, object$param$beta_kernel,object$param$tau,threshold)
    }
  }
   return(as.numeric(pred))
}

#' Predict method for the LSBP
#' 
#' 
#' Predict method for the LSBP estimated using the Gibbs sampling algorithm.
#' 
#' @param object An object of class \verb{LSBP_Gibbs}.
#' @param type String indicating the type of prediction: \verb{type="mean"},\verb{type="predictive"},\verb{type="variance"} or \verb{type="cdf"}. See details.
#' @param newdata A new data frame containing the same variables declared in \verb{Formula}. If missing, the dataset provided for estimation is used.
#' @param threshold Only needed if \verb{type="cdf"} is selected. See details.
#' @param ... Further arguments passed to or from other methods.
#' @details The method \verb{predict.LSBP_Gibbs} produces a sample of predicted values, obtained by evaluating the conditional mean of the LSBP model or the predictive distribution, using the observations contained in the \verb{newdata} data frame. 
#' 
#' If \verb{type="mean"} a sample from the posterior distribution of the LSBP mean is returned. If \verb{type="predictive"} is selected, then a sample from the predictive distribution is returned.  If \verb{type="variance"} a sample from the posterior distribution of the LSBP variance is returned.  If \verb{type="cdf"} a sample from the posterior distribution of the LSBP cumulative distribution function is returned, evaluated at \verb{threshold}.
#' @export
predict.LSBP_Gibbs <- function(object, type = "mean", newdata = NULL,threshold=NULL,...) {
  if (!(type %in% c("mean", "predictive","variance","cdf"))) 
    stop("Please provide a valid predictive type")
  
   if (is.null(newdata)) {
      X1 <- object$data$X1
      X2 <- object$data$X2
   } else {
      f <- object$call
      X1 <- model.matrix(f, data = newdata, rhs = 1)
      if (length(f)[2] == 1) {
         X2 <- X1
      } else {
         X2 <- model.matrix(f, data = newdata, rhs = 2)
      }
   }
   
   n <- NROW(X2); p_kernel <- NCOL(X1); p_mixing <- NCOL(X2)
   
   if (type == "mean") {
      pred_mean <- matrix(0, object$control$R, n)
      if (NCOL(X1) > 1) {
         # Multivariate case
         for (r in 1:object$control$R) {
            pred_mean[r, ] <- pred_mean_multi(X1, X2, matrix(object$param$beta_mixing[r, , ],object$H-1, p_mixing), 
              object$param$beta_kernel[r, , ])
         }
      } else {
         # Univariate case
         for (r in 1:object$control$R) {
            pred_mean[r, ] <- pred_mean(X2, matrix(object$param$beta_mixing[r, , ], object$H-1, p_mixing), object$param$beta_kernel[r, 
              ])
         }
      }
      return(pred_mean)
   }
   
   if (type == "predictive") {
      predictive <- matrix(0, object$control$R, n)
      if (NCOL(X1) > 1) {
         # Multivariate case
         for (r in 1:object$control$R) {
            predictive[r, ] <- predictive_multi(X1, X2, matrix(object$param$beta_mixing[r, , ],object$H-1,  p_mixing), 
              object$param$beta_kernel[r, , ], object$param$tau[r, ])
         }
      } else {
         # Univariate case
         for (r in 1:object$control$R) {
            predictive[r, ] <- predictive(X2, matrix(object$param$beta_mixing[r, , ],object$H-1,  p_mixing), object$param$beta_kernel[r, 
              ], object$param$tau[r, ])
         }
      }
      return(predictive)
   }
   
   if (type == "variance") {
     variance <- matrix(0, object$control$R, n)
     if (NCOL(X1) > 1) {
       # Multivariate case
       for (r in 1:object$control$R) {
         variance[r, ] <- pred_var_multi(X1, X2, matrix(object$param$beta_mixing[r, , ],object$H-1,  p_mixing), object$param$beta_kernel[r, , ], object$param$tau[r, ])
       }
     } else {
       # Univariate case
       for (r in 1:object$control$R) {
         variance[r, ] <- pred_var(X2, matrix(object$param$beta_mixing[r, , ],object$H-1,  p_mixing), object$param$beta_kernel[r,], object$param$tau[r, ])
       }
     }
     return(variance)
   }
   
   if (type == "cdf") {
     cdf <- matrix(0, object$control$R, n)
     if (NCOL(X1) > 1) {
       # Multivariate case
       for (r in 1:object$control$R) {
         cdf[r, ] <- pred_cdf_multi(X1, X2, matrix(object$param$beta_mixing[r, , ],object$H-1,  p_mixing), object$param$beta_kernel[r, , ], object$param$tau[r, ],threshold)
       }
     } else {
       # Univariate case
       for (r in 1:object$control$R) {
         cdf[r, ] <- pred_cdf(X2, matrix(object$param$beta_mixing[r, , ],object$H-1,  p_mixing), object$param$beta_kernel[r,], object$param$tau[r, ],threshold)
       }
     }
     return(cdf)
   }
   
}

#' Predict method for the LSBP
#' 
#' 
#' Predict method for the LSBP estimated using the Variational Bayes algorithm.
#' 
#' @param object An object of class \verb{LSBP_VB}.
#' @param type String indicating the type of prediction: \verb{type="mean"},\verb{type="predictive"},  \verb{type="variance"} or \verb{type="cdf"}. See details.
#' @param R An integer indicating the number of replications for the returned sample.
#' @param newdata A new data frame containing the same variables declared in \verb{Formula}. If missing, the dataset provided for estimation is used.
#' @param threshold Only needed if \verb{type="cdf"} is selected. See details.
#' @param ... Further arguments passed to or from other methods.
#' @details The method \verb{predict.LSBP_VB} produces a sample of predicted values, obtained by evaluating the conditional mean of the LSBP model or the predictive distribution, using the observations contained in the \verb{newdata} data frame. 
#' 
#' If \verb{type="mean"} a sample from the posterior distribution of the LSBP mean is returned. If \verb{type="predictive"} is selected, then a sample from the predictive distribution is returned.  If \verb{type="variance"} a sample from the posterior distribution of the LSBP variance is returned.  If \verb{type="cdf"} a sample from the posterior distribution of the LSBP cumulative distribution function is returned, evaluated at \verb{threshold}.
#' @export
predict.LSBP_VB <- function(object, type = "mean", R = 5000, newdata = NULL, threshold=NULL, ...) {
  if (!(type %in% c("mean", "predictive","variance","cdf"))) 
    stop("Please provide a valid predictive type")
   if (is.null(newdata)) {
      X1 <- object$data$X1
      X2 <- object$data$X2
   } else {
      f <- object$call
      X1 <- model.matrix(f, data = newdata, rhs = 1)
      if (length(f)[2] == 1) {
         X2 <- X1
      } else {
         X2 <- model.matrix(f, data = newdata, rhs = 2)
      }
   }
   
   n <- NROW(X2)
   p_kernel <- NCOL(X1)
   p_mixing <- NCOL(X2)
   beta_mixing <- array(0, c(R, object$H - 1, p_mixing))
   tau <- matrix(0, R, object$H)
   
   # Generating the parameters
   
   if (NCOL(X1) > 1) {
      beta_kernel <- array(0, c(R, object$H, p_kernel))
      for (h in 1:object$H) {
         if (h < object$H) {
            eig <- eigen(object$param$Sigma_mixing[h, , ], symmetric = TRUE)
            A1 <- t(eig$vectors) * sqrt(eig$values)
            beta_mixing[, h, ] <- t(object$param$mu_mixing[h, ] + t(matrix(rnorm(R * p_mixing), R, p_mixing) %*% 
              A1))
         }
         
         eig <- eigen(object$param$Sigma_kernel[h, , ], symmetric = TRUE)
         A1 <- t(eig$vectors) * sqrt(eig$values)
         beta_kernel[, h, ] <- t(object$param$mu_kernel[h, ] + t(matrix(rnorm(R * p_kernel), R, p_kernel) %*% 
            A1))
         tau[, h] <- rgamma(R, object$param$a_tilde[h], object$param$b_tilde[h])
      }
   } else {
      beta_kernel <- matrix(0, R, object$H)
      for (h in 1:object$H) {
         if (h < object$H) {
            eig <- eigen(object$param$Sigma_mixing[h, , ], symmetric = TRUE)
            A1 <- t(eig$vectors) * sqrt(eig$values)
            beta_mixing[, h, ] <- t(object$param$mu_mixing[h, ] + t(matrix(rnorm(R * p_mixing), R, p_mixing) %*% 
              A1))
         }
         beta_kernel[, h] <- rnorm(R, object$param$mu_kernel[h], sqrt(object$param$Sigma_kernel[h]))
         tau[, h] <- rgamma(R, object$param$a_tilde[h], object$param$b_tilde[h])
      }
      
   }
   
   if(type == "mean") {
      pred <- matrix(0, R, n)
      if (NCOL(X1) > 1) {
         for (r in 1:R) {
            pred[r, ] <- pred_mean_multi(X1, X2, matrix(beta_mixing[r, , ], object$H - 1, p_mixing), beta_kernel[r, 
              , ])
         }
      } else {
         for (r in 1:R) {
            pred[r, ] <- pred_mean(X2, matrix(beta_mixing[r, , ], object$H - 1, p_mixing), beta_kernel[r, 
              ])
         }
      }
   }
   
   if(type == "predictive") {
      pred <- matrix(0, R, n)
      if (NCOL(X1) > 1) {
         for (r in 1:R) {
            pred[r, ] <- predictive_multi(X1, X2, matrix(beta_mixing[r, , ], object$H - 1, p_mixing), 
              beta_kernel[r, , ], tau[r, ])
         }
      } else {
         for (r in 1:R) {
            pred[r, ] <- predictive(X2, matrix(beta_mixing[r, , ], object$H - 1, p_mixing), beta_kernel[r, 
              ], tau[r, ])
         }
      }
   }
   
   if(type == "variance") {
     pred <- matrix(0, R, n)
     if (NCOL(X1) > 1) {
       for (r in 1:R) {
         pred[r, ] <- pred_var_multi(X1, X2, matrix(beta_mixing[r, , ], object$H - 1, p_mixing), 
                                       beta_kernel[r, , ], tau[r, ])
       }
     } else {
       for (r in 1:R) {
         pred[r, ] <- pred_var(X2, matrix(beta_mixing[r, , ], object$H - 1, p_mixing), beta_kernel[r, 
                                                                                                     ], tau[r, ])
       }
     }
   }
   
   if(type == "cdf") {
     pred <- matrix(0, R, n)
     if (NCOL(X1) > 1) {
       for (r in 1:R) {
         pred[r, ] <- pred_cdf_multi(X1, X2, matrix(beta_mixing[r, , ], object$H - 1, p_mixing), 
                                     beta_kernel[r, , ], tau[r, ],threshold)
       }
     } else {
       for (r in 1:R) {
         pred[r, ] <- pred_cdf(X2, matrix(beta_mixing[r, , ], object$H - 1, p_mixing), beta_kernel[r, 
                                                                                                   ], tau[r, ], threshold)
       }
     }
   }
   
   return(pred)
}

#' Conditional density for a LSBP model
#' 
#' 
#' Evaluate the conditional density of a LSBP given the parameters.
#' 
#' @param y A value of which the conditional density has to be computed
#' @param X1 A \verb{n x p_kernel} design matrix for the kernel
#' @param X2 A \verb{n x p_mixing} design matrix for the stick-breaking weights
#' @param beta A \verb{H x p_kernel} dimensional matrix of coefficients for the linear predictor of the kernel
#' @param alpha A \verb{H-1 x p_mixing} dimensional matrix of coefficients for the linear predictor of the stick-breaking weights
#' @param tau A \verb{H} dimensional vector of coefficients for the kernel precision
#' @details The function \verb{LSBP_density} evaluates the conditional density of \verb{y} given the parameters
#' @export
LSBP_density <- function(y, X1, X2, beta, alpha, tau){
  LSBP_density_C(y, X1, X2, beta, alpha, tau)
}