Lint using styler

nathansam · May 9, 2022 · d0cc4e8 · d0cc4e8
1 parent 89a89bc
commit d0cc4e8
Show file tree

Hide file tree

Showing 18 changed files with 3,114 additions and 2,556 deletions.
diff --git a/R/Internal_Codes.R b/R/Internal_Codes.R
diff --git a/R/LogExponentialPower.R b/R/LogExponentialPower.R
diff --git a/R/LogLaplace.R b/R/LogLaplace.R
@@ -1,5 +1,5 @@
 ################################################################################
-#####################IMPLEMENTATION OF THE LOG-LAPLACE MODEL####################
+##################### IMPLEMENTATION OF THE LOG-LAPLACE MODEL####################
 ################################################################################
 #' @title MCMC algorithm for the log-Laplace model
 #' @description  Adaptive Metropolis-within-Gibbs algorithm with univariate
@@ -20,8 +20,10 @@
 #' # This is only an illustration. Run longer chains for more accurate
 #' # estimations.
 #'
-#' LLAP <- MCMC_LLAP(N = 1000, thin = 20, burn = 40, Time = cancer[, 1],
-#'                   Cens = cancer[, 2], X = cancer[, 3:11])
+#' LLAP <- MCMC_LLAP(
+#'   N = 1000, thin = 20, burn = 40, Time = cancer[, 1],
+#'   Cens = cancer[, 2], X = cancer[, 3:11]
+#' )
 #'
 #' @export
 MCMC_LLAP <- function(N,
@@ -42,18 +44,20 @@ MCMC_LLAP <- function(N,
   if (is.null(beta0)) beta0 <- beta.sample(n = ncol(X))
   if (is.null(sigma20)) sigma20 <- sigma2.sample()
 
-  MCMC.param.check(N,
-                   thin,
-                   burn,
-                   Time,
-                   Cens,
-                   X,
-                   beta0,
-                   sigma20,
-                   prior,
-                   set,
-                   eps_l,
-                   eps_r)
+  MCMC.param.check(
+    N,
+    thin,
+    burn,
+    Time,
+    Cens,
+    X,
+    beta0,
+    sigma20,
+    prior,
+    set,
+    eps_l,
+    eps_r
+  )
 
   k <- length(beta0)
   n <- length(Time)
@@ -72,7 +76,7 @@ MCMC_LLAP <- function(N,
   logt <- matrix(rep(0, times = (N.aux + 1) * n), ncol = n)
   logt[1, ] <- log(Time)
   lambda <- matrix(rep(0, times = (N.aux + 1) * n), ncol = n)
-  lambda[1, ] <- (stats::rgamma(n, shape = 2, rate = 2)) ^ (-1)
+  lambda[1, ] <- (stats::rgamma(n, shape = 2, rate = 2))^(-1)
 
   beta.aux <- beta[1, ]
   sigma2.aux <- sigma2[1]
@@ -97,8 +101,9 @@ MCMC_LLAP <- function(N,
       (logt.aux - X %*% beta.aux)
     if (rate.aux > 0 & is.na(rate.aux) == FALSE) {
       sigma2.aux <- (stats::rgamma(1,
-                                   shape = shape.aux,
-                                   rate = rate.aux)) ^ (-1)
+        shape = shape.aux,
+        rate = rate.aux
+      ))^(-1)
       if (sigma2.aux < 0) {
         cat("Negative sigma2")
       }
@@ -107,9 +112,11 @@ MCMC_LLAP <- function(N,
     if ((iter - 1) %% Q == 0 && iter - 1 <= burn) {
       mu.aux <- sqrt(sigma2.aux) / abs(logt.aux - X %*% beta.aux)
       if (sum(is.na(mu.aux)) == 0) {
-        draw.aux <- VGAM::rinv.gaussian(n = n,
-                                        mu = mu.aux,
-                                        lambda = rep(1, times = n))
+        draw.aux <- VGAM::rinv.gaussian(
+          n = n,
+          mu = mu.aux,
+          lambda = rep(1, times = n)
+        )
         lambda.aux <- I(draw.aux > 0) * draw.aux +
           (1 - I(draw.aux > 0)) * lambda.aux
         lambda.aux.change <- TRUE
@@ -118,14 +125,16 @@ MCMC_LLAP <- function(N,
       }
     }
 
-    logt.aux <- logt.update.SMLN(Time,
-                                 Cens,
-                                 X,
-                                 beta.aux,
-                                 sigma2.aux / lambda.aux,
-                                 set,
-                                 eps_l,
-                                 eps_r)
+    logt.aux <- logt.update.SMLN(
+      Time,
+      Cens,
+      X,
+      beta.aux,
+      sigma2.aux / lambda.aux,
+      set,
+      eps_l,
+      eps_r
+    )
 
     if (iter %% thin == 0) {
       beta[iter / thin + 1, ] <- beta.aux
@@ -150,7 +159,7 @@ MCMC_LLAP <- function(N,
 
   if (burn > 0) {
     burn.period <- 1:(burn / thin)
-    chain <- chain [-burn.period, ]
+    chain <- chain[-burn.period, ]
   }
 
   return(chain)
@@ -167,8 +176,10 @@ MCMC_LLAP <- function(N,
 #' # This is only an illustration. Run longer chains for more accurate
 #' # estimations.
 #'
-#' LLAP <- MCMC_LLAP(N = 1000, thin = 20, burn = 40, Time = cancer[, 1],
-#'                   Cens = cancer[, 2], X = cancer[, 3:11])
+#' LLAP <- MCMC_LLAP(
+#'   N = 1000, thin = 20, burn = 40, Time = cancer[, 1],
+#'   Cens = cancer[, 2], X = cancer[, 3:11]
+#' )
 #'
 #' @export
 LML_LLAP <- function(thin,
@@ -201,53 +212,62 @@ LML_LLAP <- function(thin,
 
   # LIKELIHOOD ORDINATE
   LL.ord <- log.lik.LLAP(Time,
-                         Cens,
-                         X,
-                         beta = beta.star,
-                         sigma2 = sigma2.star,
-                         set,
-                         eps_l,
-                         eps_r)
+    Cens,
+    X,
+    beta = beta.star,
+    sigma2 = sigma2.star,
+    set,
+    eps_l,
+    eps_r
+  )
   cat("Likelihood ordinate ready!\n")
 
   # PRIOR ORDINATE
-  LP.ord <- prior_LN(beta = beta.star,
-                     sigma2 = sigma2.star,
-                     prior = prior,
-                     logs = TRUE)
+  LP.ord <- prior_LN(
+    beta = beta.star,
+    sigma2 = sigma2.star,
+    prior = prior,
+    logs = TRUE
+  )
   cat("Prior ordinate ready!\n")
 
   # POSTERIOR ORDINATE - sigma 2
   shape <- (n + 2 * p - 2) / 2
   po.sigma2 <- rep(0, times = N)
   for (i in 1:N) {
-    aux1 <- as.vector(t(chain[i,
-                              (k + 2 + n):(k + 1 + 2 * n)])) - X %*%
+    aux1 <- as.vector(t(chain[
+      i,
+      (k + 2 + n):(k + 1 + 2 * n)
+    ])) - X %*%
       as.vector(chain[i, 1:k])
     aux2 <- diag(as.vector(t(chain[i, (k + 2):(k + 1 + n)])))
     rate.aux <- as.numeric(0.5 * t(aux1) %*% aux2 %*% aux1)
-    po.sigma2[i] <- MCMCpack::dinvgamma(sigma2.star, shape = shape,
-                                        scale = rate.aux)
+    po.sigma2[i] <- MCMCpack::dinvgamma(sigma2.star,
+      shape = shape,
+      scale = rate.aux
+    )
   }
   PO.sigma2 <- mean(po.sigma2)
   cat("Posterior ordinate sigma2 ready!\n")
 
   # POSTERIOR ORDINATE - beta
   logt0 <- t(chain[N, (n + k + 2):(2 * n + k + 1)])
-  chain.beta <- MCMCR.sigma2.LLAP(N = N * thin,
-                                  thin = thin,
-                                  Q,
-                                  Time,
-                                  Cens,
-                                  X,
-                                  beta0 = t(chain[N, 1:k]),
-                                  sigma20 = sigma2.star,
-                                  logt0 = logt0,
-                                  lambda0 = t(chain[N, (k + 2):(n + k + 1)]),
-                                  prior = prior,
-                                  set,
-                                  eps_l,
-                                  eps_r)
+  chain.beta <- MCMCR.sigma2.LLAP(
+    N = N * thin,
+    thin = thin,
+    Q,
+    Time,
+    Cens,
+    X,
+    beta0 = t(chain[N, 1:k]),
+    sigma20 = sigma2.star,
+    logt0 = logt0,
+    lambda0 = t(chain[N, (k + 2):(n + k + 1)]),
+    prior = prior,
+    set,
+    eps_l,
+    eps_r
+  )
   cat("Reduced chain.beta ready!\n")
 
   po.beta <- rep(0, times = N)
@@ -256,10 +276,11 @@ LML_LLAP <- function(thin,
     aux1.beta <- solve(t(X) %*% aux0.beta %*% X)
     aux2.beta <- aux1.beta %*% t(X) %*% aux0.beta
     mu.aux <- as.vector(aux2.beta %*%
-                          ((chain.beta[(i + 1), (n + k + 1):(2 * n + k)])))
+      ((chain.beta[(i + 1), (n + k + 1):(2 * n + k)])))
     po.beta[i] <- mvtnorm::dmvnorm(beta.star,
-                                   mean = mu.aux,
-                                   sigma = sigma2.star * aux1.beta)
+      mean = mu.aux,
+      sigma = sigma2.star * aux1.beta
+    )
   }
   PO.beta <- mean(po.beta)
   cat("Posterior ordinate beta ready!\n")
@@ -271,11 +292,13 @@ LML_LLAP <- function(thin,
   # MARGINAL LOG-LIKELIHOOD
   LML <- LL.ord + LP.ord - LPO.sigma2 - LPO.beta
 
-  return(list(LL.ord = LL.ord,
-              LP.ord = LP.ord,
-              LPO.sigma2 = LPO.sigma2,
-              LPO.beta = LPO.beta,
-              LML = LML))
+  return(list(
+    LL.ord = LL.ord,
+    LP.ord = LP.ord,
+    LPO.sigma2 = LPO.sigma2,
+    LPO.beta = LPO.beta,
+    LML = LML
+  ))
 }
 
 #' @title Deviance information criterion for the log-Laplace model
@@ -291,10 +314,14 @@ LML_LLAP <- function(thin,
 #' # This is only an illustration. Run longer chains for more accurate
 #' # estimations.
 #'
-#' LLAP <- MCMC_LLAP(N = 1000, thin = 20, burn = 40, Time = cancer[, 1],
-#'                   Cens = cancer[, 2], X = cancer[, 3:11])
-#' LLAP.DIC <- DIC_LLAP(Time = cancer[, 1], Cens = cancer[, 2],
-#'                      X = cancer[, 3:11], chain = LLAP)
+#' LLAP <- MCMC_LLAP(
+#'   N = 1000, thin = 20, burn = 40, Time = cancer[, 1],
+#'   Cens = cancer[, 2], X = cancer[, 3:11]
+#' )
+#' LLAP.DIC <- DIC_LLAP(
+#'   Time = cancer[, 1], Cens = cancer[, 2],
+#'   X = cancer[, 3:11], chain = LLAP
+#' )
 #'
 #' @export
 DIC_LLAP <- function(Time,
@@ -311,23 +338,25 @@ DIC_LLAP <- function(Time,
 
   for (iter in 1:N) {
     LL[iter] <- log.lik.LLAP(Time,
-                             Cens, X,
-                             beta = as.vector(chain[iter, 1:k]),
-                             sigma2 = chain[iter, k + 1],
-                             set = set,
-                             eps_l,
-                             eps_r)
+      Cens, X,
+      beta = as.vector(chain[iter, 1:k]),
+      sigma2 = chain[iter, k + 1],
+      set = set,
+      eps_l,
+      eps_r
+    )
   }
 
   aux <- apply(chain[, 1:(k + 1)], 2, "median")
   pd <- -2 * mean(LL) + 2 * log.lik.LLAP(Time,
-                                         Cens,
-                                         X,
-                                         beta = aux[1:k],
-                                         sigma2 = aux[k + 1],
-                                         set = set,
-                                         eps_l,
-                                         eps_r)
+    Cens,
+    X,
+    beta = aux[1:k],
+    sigma2 = aux[k + 1],
+    set = set,
+    eps_l,
+    eps_r
+  )
   pd.aux <- k + 1
 
   DIC <- -2 * mean(LL) + pd
@@ -354,10 +383,14 @@ DIC_LLAP <- function(Time,
 #' # This is only an illustration. Run longer chains for more accurate
 #' # estimations.
 #'
-#' LLAP <- MCMC_LLAP(N = 1000, thin = 20, burn = 40, Time = cancer[, 1],
-#'                   Cens = cancer[, 2], X = cancer[, 3:11])
-#' LLAP.CD <- CaseDeletion_LLAP(Time = cancer[, 1], Cens = cancer[, 2],
-#'                              X = cancer[, 3:11], chain = LLAP)
+#' LLAP <- MCMC_LLAP(
+#'   N = 1000, thin = 20, burn = 40, Time = cancer[, 1],
+#'   Cens = cancer[, 2], X = cancer[, 3:11]
+#' )
+#' LLAP.CD <- CaseDeletion_LLAP(
+#'   Time = cancer[, 1], Cens = cancer[, 2],
+#'   X = cancer[, 3:11], chain = LLAP
+#' )
 #'
 #' @export
 CaseDeletion_LLAP <- function(Time,
@@ -379,13 +412,14 @@ CaseDeletion_LLAP <- function(Time,
     aux2 <- rep(0, times = N)
     for (ITER in 1:N) {
       aux2[ITER] <- log.lik.LLAP(Time[s],
-                                 Cens[s],
-                                 X[s, ],
-                                 beta = as.vector(chain[ITER, 1:k]),
-                                 sigma2 = chain[ITER, (k + 1)],
-                                 set,
-                                 eps_l,
-                                 eps_r)
+        Cens[s],
+        X[s, ],
+        beta = as.vector(chain[ITER, 1:k]),
+        sigma2 = chain[ITER, (k + 1)],
+        set,
+        eps_l,
+        eps_r
+      )
       aux1[ITER] <- exp(-aux2[ITER])
     }
     logCPO[s] <- -log(mean(aux1))
@@ -420,9 +454,11 @@ CaseDeletion_LLAP <- function(Time,
 #' # This is only an illustration. Run longer chains for more accurate
 #' # estimations.
 #'
-#' LLAP <- MCMC_LLAP(N = 1000, thin = 20, burn = 40, Time = cancer[, 1],
-#'                   Cens = cancer[, 2], X = cancer[, 3:11])
-#' LLAP.outlier <- BF_lambda_obs_LLAP(1,1, X = cancer[, 3:11], chain = LLAP)
+#' LLAP <- MCMC_LLAP(
+#'   N = 1000, thin = 20, burn = 40, Time = cancer[, 1],
+#'   Cens = cancer[, 2], X = cancer[, 3:11]
+#' )
+#' LLAP.outlier <- BF_lambda_obs_LLAP(1, 1, X = cancer[, 3:11], chain = LLAP)
 #'
 #' @export
 BF_lambda_obs_LLAP <- function(obs, ref, X, chain) {
@@ -435,13 +471,15 @@ BF_lambda_obs_LLAP <- function(obs, ref, X, chain) {
 
   for (j in 1:N) {
     aux1[j] <- stats::dnorm(chain[j, (obs + n + k + 1)],
-                            mean = (X[obs, ]) %*%
-                              as.vector(chain[j, 1:k]),
-                            sd = sqrt(chain[j, (k + 1)] / ref))
+      mean = (X[obs, ]) %*%
+        as.vector(chain[j, 1:k]),
+      sd = sqrt(chain[j, (k + 1)] / ref)
+    )
     aux2[j] <- VGAM::dlaplace(chain[j, (obs + n + k + 1)],
-                              location = (X[obs, ]) %*%
-                                as.vector(chain[j, 1:k]),
-                              scale = sqrt(chain[j, (k + 1)]))
+      location = (X[obs, ]) %*%
+        as.vector(chain[j, 1:k]),
+      scale = sqrt(chain[j, (k + 1)])
+    )
   }
 
   aux <- mean(aux1 / aux2)