Skip to content

Commit

Permalink
first test of ZIPLN on real data
Browse files Browse the repository at this point in the history
  • Loading branch information
@jchiquet
jchiquet committed Jan 18, 2024
1 parent 706120d commit b971360
Show file tree
Hide file tree
Showing 2 changed files with 94 additions and 12 deletions.
26 changes: 26 additions & 0 deletions inst/case_studies/oaks_tree.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,32 @@ rbind(
as.data.frame(row.names = c("full", "diagonal", "spherical")) %>%
knitr::kable()

## ZIPLN
system.time(myZIPLN_full_covar <- ZIPLN(Abundance ~ 0 + tree + offset(log(Offset)) | 0 + tree, data = oaks, control = ZIPLN_param(covariance = "full")))
system.time(myZIPLN_full_row <- ZIPLN(Abundance ~ 0 + tree + offset(log(Offset)), zi = "row", data = oaks, control = ZIPLN_param(covariance = "full")))
system.time(myZIPLN_full_col <- ZIPLN(Abundance ~ 0 + tree + offset(log(Offset)), zi = "col", data = oaks, control = ZIPLN_param(covariance = "full")))
system.time(myZIPLN_full_single <- ZIPLN(Abundance ~ 0 + tree + offset(log(Offset)), data = oaks, control = ZIPLN_param(covariance = "full")))
system.time(myZIPLN_diagonal_covar <- ZIPLN(Abundance ~ 0 + tree + offset(log(Offset)) | 0 + tree, data = oaks, control = ZIPLN_param(covariance = "diagonal")))
system.time(myZIPLN_diagonal_row <- ZIPLN(Abundance ~ 0 + tree + offset(log(Offset)), zi = "row", data = oaks, control = ZIPLN_param(covariance = "diagonal")))
system.time(myZIPLN_diagonal_col <- ZIPLN(Abundance ~ 0 + tree + offset(log(Offset)), zi = "col", data = oaks, control = ZIPLN_param(covariance = "diagonal")))
system.time(myZIPLN_diagonal_single <- ZIPLN(Abundance ~ 0 + tree + offset(log(Offset)), data = oaks, control = ZIPLN_param(covariance = "diagonal")))

rbind(
myZIPLN_full_single$criteria,
myZIPLN_full_row$criteria,
myZIPLN_full_col$criteria,
myZIPLN_full_covar$criteria,
myZIPLN_diagonal_single$criteria,
myZIPLN_diagonal_row$criteria,
myZIPLN_diagonal_col$criteria,
myZIPLN_diagonal_covar$criteria
) %>%
as.data.frame(row.names = c("ZIPLN full single", "ZIPLN full column prob", "ZIPLN full row prob", "ZIPLN full covar prob",
"ZIPLN diagonal single", "ZIPLN diagonal column prob", "ZIPLN diagonal row prob", "ZIPLN diagonal covar prob")) %>%
knitr::kable()



## Discriminant Analysis with LDA
myLDA_tree <- PLNLDA(Abundance ~ 1 + offset(log(Offset)), grouping = tree, data = oaks)
plot(myLDA_tree)
Expand Down
80 changes: 68 additions & 12 deletions inst/simus_ZIPLN/essai_ZIPLN.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,9 +38,9 @@ Omega_star <- solve(Sigma_star)
## PLN part: one intercept/grand mean with value equal to 1
B_star <- matrix(2, 1, p)
## ZI part: proba of 0.1 for being a 0 from the ZI component
B0_star <- matrix(logit(0.2), 1, p)
B0_star <- matrix(logit(0.5), 1, p)

vec_n <- c(25, 50, 75, 100)
vec_n <- c(50, 100, 200)
one_simu <- function(i) {
cat(i)
n <- max(vec_n)*2
Expand All @@ -51,24 +51,80 @@ one_simu <- function(i) {
Y <- rZIPLN(n, mu = mu_star, Sigma = Sigma_star, Pi = Pi_star)
Y <- Y[rowSums(Y) > 0, ]
X <- X[rowSums(Y) > 0]
err_ZIPLN <- sapply(vec_n, function(n_) {

err_ZIPLN_simple <- sapply(vec_n, function(n_) {
Y_ <- Y[1:n_, ]; X_ <- X[1:n_]
myZIPLN <- tryCatch(
ZIPLN(Y_ ~ 0 + X_ + offset(log(rowSums(Y_))), control = ZIPLN_param(trace = 0)),
error = function(e) {e})
if (is(myZIPLN, "error")) {
res <- c(pred_Y = NA, rmse_B = NA, rmse_Omega = NA)
} else {
res <- c(pred_Y = sqrt(mean((myZIPLN$fitted - Y_)^2)),
rmse_B = sqrt(mean((myZIPLN$model_par$B - B_star)^2)),
rmse_Omega = sqrt(mean((myZIPLN$model_par$Omega - Omega_star)^2)))
}
res
})

err_ZIPLN_row <- sapply(vec_n, function(n_) {
Y_ <- Y[1:n_, ]; X_ <- X[1:n_]
myZIPLN <- ZIPLN(Y_ ~ 0 + X_ + offset(log(rowSums(Y_))) | X_ , control = ZIPLN_param(trace = 0))
c(pred_Y = sqrt(mean((myZIPLN$fitted - Y_)^2)),
rmse_B = sqrt(mean((myZIPLN$model_par$B[2, ] - B_star)^2)),
myZIPLN <- tryCatch(
ZIPLN(Y_ ~ 0 + X_ + offset(log(rowSums(Y_))), zi = "row", control = ZIPLN_param(trace = 0)),
error = function(e) {e})
if (is(myZIPLN, "error")) {
res <- c(pred_Y = NA, rmse_B = NA, rmse_Omega = NA)
} else {
res <- c(pred_Y = sqrt(mean((myZIPLN$fitted - Y_)^2)),
rmse_B = sqrt(mean((myZIPLN$model_par$B - B_star)^2)),
rmse_Omega = sqrt(mean((myZIPLN$model_par$Omega - Omega_star)^2)))
}
res
})

err_ZIPLN_col <- sapply(vec_n, function(n_) {
Y_ <- Y[1:n_, ]; X_ <- X[1:n_]
myZIPLN <- tryCatch(
myZIPLN <- ZIPLN(Y_ ~ 0 + X_ + offset(log(rowSums(Y_))), zi = "col", control = ZIPLN_param(trace = 0)),
error = function(e) {e})
if (is(myZIPLN, "error")) {
res <- c(pred_Y = NA, rmse_B = NA, rmse_Omega = NA)
} else {
res <- c(pred_Y = sqrt(mean((myZIPLN$fitted - Y_)^2)),
rmse_B = sqrt(mean((myZIPLN$model_par$B - B_star)^2)),
rmse_Omega = sqrt(mean((myZIPLN$model_par$Omega - Omega_star)^2)))
}
res
})

err_ZIPLN_covar <- sapply(vec_n, function(n_) {
Y_ <- Y[1:n_, ]; X_ <- X[1:n_]

myZIPLN <- tryCatch(
myZIPLN <- ZIPLN(Y_ ~ 0 + X_ + offset(log(rowSums(Y_))) | X_ , control = ZIPLN_param(trace = 0)),
error = function(e) {e})
if (is(myZIPLN, "error")) {
res <- c(pred_Y = NA, rmse_B = NA, rmse_Omega = NA)
} else {
res <- c(pred_Y = sqrt(mean((myZIPLN$fitted - Y_)^2)),
rmse_B = sqrt(mean((myZIPLN$model_par$B - B_star)^2)),
rmse_Omega = sqrt(mean((myZIPLN$model_par$Omega - Omega_star)^2)))
}
res
})

err_PLN <- sapply(vec_n, function(n_) {
Y_ <- Y[1:n_, ]; X_ <- X[1:n_]
myPLN <- PLN(Y_ ~ 1 + X_ + offset(log(rowSums(Y_))) , control = PLN_param(trace = 0))
myPLN <- PLN(Y_ ~ 0 + X_ + offset(log(rowSums(Y_))) , control = PLN_param(trace = 0))
c(pred = sqrt(mean((myPLN$fitted - Y_)^2)),
rmse_B = sqrt(mean((myPLN$model_par$B[2, ] - B_star)^2)),
rmse_B = sqrt(mean((myPLN$model_par$B - B_star)^2)),
rmse_Sigma = sqrt(mean((myPLN$model_par$Omega - Omega_star)^2))
)
})
data.frame(rbind(t(err_ZIPLN), t(err_PLN)),
method = rep(c("ZIPLN","PLN"), each = length(vec_n)),
n = c(vec_n, vec_n), simu = rep(i, 2*length(vec_n)))

data.frame(rbind(t(err_ZIPLN_simple), t(err_ZIPLN_row), t(err_ZIPLN_col), t(err_ZIPLN_covar), t(err_PLN)),
method = rep(c("ZIPLN_simple", "ZIPLN_row", "ZIPLN_col", "ZIPLN_covar", "PLN"), each = length(vec_n)),
n = rep(vec_n, 5), simu = rep(i, 5*length(vec_n)))
}

res <- do.call(rbind, lapply(1:50, one_simu))
Expand All @@ -77,7 +133,7 @@ p <- ggplot(res) + aes(x = factor(n), y = pred_Y, fill = factor(method)) + geom_
scale_y_log10() + ylim(c(0,2))
p

p <- ggplot(res) + aes(x = factor(n), y = rmse_B, fill = factor(method)) + geom_violin() + theme_bw() + scale_y_log10() + ylim(c(0,3))
p <- ggplot(res) + aes(x = factor(n), y = rmse_B, fill = factor(method)) + geom_violin() + theme_bw() + scale_y_log10() + ylim(c(2.75,3))
p

p <- ggplot(res) + aes(x = factor(n), y = rmse_Omega, fill = factor(method)) + geom_violin() + theme_bw() + scale_y_log10() + ylim(c(0,0.5))
Expand Down

0 comments on commit b971360

Please sign in to comment.