Skip to content

Commit

Permalink
Change initialization in variance_*() to mimic a brand new call to …
Browse files Browse the repository at this point in the history 
…`PLN()`, based on problems reported in #112.
  • Loading branch information
@mahendra-mariadassou
mahendra-mariadassou committed Nov 17, 2023
1 parent df51164 commit c740f2d
Showing 1 changed file with 28 additions and 18 deletions.
46 changes: 28 additions & 18 deletions R/PLNfit-class.R
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,7 @@
## %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
## CLASS PLNfit #######################################
## %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

#' An R6 Class to represent a PLNfit in a standard, general framework
#'
#' @description The function [PLN()] fit a model which is an instance of a object with class [`PLNfit`].
Expand Down Expand Up @@ -207,15 +211,15 @@ PLNfit <- R6Class(
variance_jackknife = function(Y, X, O, w, config = config_default_nlopt) {
cat("\n Computing jackknife variance estimator", sep = "\n")
jacks <- future.apply::future_lapply(seq_len(self$n), function(i) {
# cat(paste0("Jackknife estimate ", i, "/", self$n), sep = "\n")
data <- list(Y = Y[-i, , drop = FALSE],
X = X[-i, , drop = FALSE],
O = O[-i, , drop = FALSE],
w = w[-i])
args <- list(data = data,
params = list(B = private$B * (1 + matrix(runif(min = -0.5, max = 0.5, n = self$d * self$p), self$d, self$p)),
M = matrix(0, self$n-1, self$p),
S = private$S[-i, , drop = FALSE]),
# params = list(B = private$B,
# M = matrix(0, self$n-1, self$p),
# S = private$S[-i, , drop = FALSE]),
params = do.call(starting_point, data),
config = config)
optim_out <- do.call(private$optimizer$main, args)
optim_out[c("B", "Omega")]
Expand All @@ -227,7 +231,6 @@ PLNfit <- R6Class(
B_hat <- private$B[,] ## strips attributes while preserving names
attr(private$B, "bias") <- (self$n - 1) * (B_jack - B_hat)
attr(private$B, "variance_jackknife") <- (self$n - 1) / self$n * var_jack
attr(private$B, "estimates_jackknife") <- jacks

Omega_jack <- jacks %>% map("Omega") %>% reduce(`+`) / self$n
var_jack <- jacks %>% map("Omega") %>% map(~( (. - Omega_jack)^2)) %>% reduce(`+`) %>%
Expand All @@ -245,7 +248,8 @@ PLNfit <- R6Class(
O = O[resample, , drop = FALSE],
w = w[resample])
args <- list(data = data,
params = list(B = private$B, M = matrix(0,self$n,self$p), S = private$S[resample, ]),
# params = list(B = private$B, M = matrix(0,self$n,self$p), S = private$S[resample, ]),
params = do.call(starting_point, data),
config = config)
optim_out <- do.call(private$optimizer$main, args)
optim_out[c("B", "Omega", "monitoring")]
Expand Down Expand Up @@ -306,10 +310,14 @@ PLNfit <- R6Class(
private$S <- control$inception$var_par$S
} else {
if (control$trace > 1) cat("\n Use LM after log transformation to define the inceptive model")
fits <- lm.fit(weights * covariates, weights * log((1 + responses)/exp(offsets)))
private$B <- matrix(fits$coefficients, d, p)
private$M <- matrix(fits$residuals, n, p)
private$S <- matrix(.1, n, p)
# fits <- lm.fit(weights * covariates, weights * log((1 + responses)/exp(offsets)))
# private$B <- matrix(fits$coefficients, d, p)
# private$M <- matrix(fits$residuals, n, p)
# private$S <- matrix(.1, n, p)
start_point <- starting_point(Y = responses, X = covariates, O = offsets, w = weights)
private$B <- start_point$B
private$M <- start_point$M
private$S <- start_point$S
}
private$optimizer$main <- ifelse(control$backend == "nlopt", nlopt_optimize, private$torch_optimize)
private$optimizer$vestep <- nlopt_optimize_vestep
Expand Down Expand Up @@ -600,7 +608,7 @@ PLNfit <- R6Class(
)

## %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
## CLASS PLNfit_diagonal
## CLASS PLNfit_diagonal ##############################
## %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

#' An R6 Class to represent a PLNfit in a standard, general framework, with diagonal residual covariance
Expand Down Expand Up @@ -682,7 +690,7 @@ PLNfit_diagonal <- R6Class(
)

## %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
## CLASS PLNfit_spherical
## CLASS PLNfit_spherical ############################
## %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

#' An R6 Class to represent a PLNfit in a standard, general framework, with spherical residual covariance
Expand Down Expand Up @@ -762,7 +770,7 @@ PLNfit_spherical <- R6Class(
)

## %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
## CLASS PLNfit_fixedcov
## CLASS PLNfit_fixedcov #############################
## %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

#' An R6 Class to represent a PLNfit in a standard, general framework, with fixed (inverse) residual covariance
Expand Down Expand Up @@ -855,11 +863,13 @@ PLNfit_fixedcov <- R6Class(

variance_jackknife = function(Y, X, O, w, config = config_default_nlopt) {
jacks <- future.apply::future_lapply(seq_len(self$n), function(i) {
args <- list(Y = Y[-i, , drop = FALSE],
data <- list(Y = Y[-i, , drop = FALSE],
X = X[-i, , drop = FALSE],
O = O[-i, , drop = FALSE],
w = w[-i],
params = list(B = private$B, Omega = private$Omega, M = private$M[-i, ], S = private$S[-i, ]),
w = w[-i])
args <- list(data = data,
# params = list(B = private$B, Omega = private$Omega, M = private$M[-i, ], S = private$S[-i, ]),
params = do.call(starting_point, data),
config = config)
optim_out <- do.call(private$optimizer$main, args)
optim_out[c("B", "Omega")]
Expand Down Expand Up @@ -912,8 +922,8 @@ PLNfit_fixedcov <- R6Class(
)

## %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
## CLASS PLNfit_genetprior
# ## %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
## CLASS PLNfit_genetprior ############################
# ## %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#
# #' An R6 Class to represent a PLNfit in a standard, general framework, with residual covariance modelling
# #' motivatived by population genetics
Expand Down

0 comments on commit c740f2d

Please sign in to comment.