diff --git a/DESCRIPTION b/DESCRIPTION
index eed5f0a2..f53a5069 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: PLNmodels
 Title: Poisson Lognormal Models
-Version: 1.0.4-0200
+Version: 1.0.4-0300
 Authors@R: c(
   person("Julien", "Chiquet", role = c("aut", "cre"), email = "julien.chiquet@inrae.fr",
       comment = c(ORCID = "0000-0002-3629-3429")),
diff --git a/NAMESPACE b/NAMESPACE
index cd0efa86..b21b7021 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -43,6 +43,7 @@ export(PLNmixture_param)
 export(PLNnetwork)
 export(PLNnetwork_param)
 export(coefficient_path)
+export(compute_PLN_starting_point)
 export(compute_offset)
 export(extract_probs)
 export(getBestModel)
diff --git a/NEWS.md b/NEWS.md
index 6484bf95..428c9259 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -2,6 +2,7 @@
 
 * Update documentation of PLN*_param() functions to include torch optimization parameters
 * Add (somehow) explicit error message when torch convergence fails
+* Change initialization in `variance_jackknife()` and `variance_bootstrap()` to prevent estimation recycling, results from those functions are now comparable to doing jackknife / bootstrap "by hand". 
 
 # PLNmodels 1.0.4 (2023-08-24)
 
diff --git a/R/PLNfit-class.R b/R/PLNfit-class.R
index 715d0499..9b9f9e9a 100644
--- a/R/PLNfit-class.R
+++ b/R/PLNfit-class.R
@@ -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`].
@@ -211,7 +215,10 @@ PLNfit <- R6Class(
                      O = O[-i, , drop = FALSE],
                      w = w[-i])
         args <- list(data = data,
-                     params = list(B = private$B, M = matrix(0, self$n-1, self$p), S = private$S[-i, ]),
+                     # params = list(B = private$B,
+                     #               M = matrix(0, self$n-1, self$p),
+                     #               S = private$S[-i, , drop = FALSE]),
+                     params = do.call(compute_PLN_starting_point, data),
                      config = config)
         optim_out <- do.call(private$optimizer$main, args)
         optim_out[c("B", "Omega")]
@@ -240,7 +247,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(compute_PLN_starting_point, data),
                      config = config)
         optim_out <- do.call(private$optimizer$main, args)
         optim_out[c("B", "Omega", "monitoring")]
@@ -301,10 +309,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 <- compute_PLN_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
@@ -595,7 +607,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
@@ -677,7 +689,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
@@ -757,7 +769,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
@@ -850,11 +862,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(compute_PLN_starting_point, data),
                      config = config)
         optim_out <- do.call(private$optimizer$main, args)
         optim_out[c("B", "Omega")]
@@ -907,8 +921,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
diff --git a/R/utils.R b/R/utils.R
index 9bc10dfe..87e59712 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -252,3 +252,39 @@ create_parameters <- function(
        Sigma  = sigma * toeplitz(x = rho^seq(0, p-1)),
        depths = depths)
 }
+
+#' Helper function for PLN initialization.
+#'
+#' @description
+#' Barebone function to compute starting points for B, M and S when fitting a PLN. Mostly intended for internal use.
+#'
+#' @param Y Response count matrix
+#' @param X Covariate matrix
+#' @param O Offset matrix (in log-scale)
+#' @param w Weight vector (defaults to 1)
+#' @param s Scale parameter for S (defaults to 0.1)
+#' @return a named list of starting values for model parameter B and variational parameters M and S used in the iterative optimization algorithm of [PLN()]
+#'
+#' @details The default strategy to estimate B and M is to fit a linear model with covariates `X` to the response count matrix (after adding a pseudocount of 1, scaling by the offset and taking the log). The regression matrix is used to initialize `B` and the residuals to initialize `M`. `S` is initialized as a constant conformable matrix with value `s`.
+#'
+#' @rdname compute_PLN_starting_point
+#' @examples
+#' \dontrun{
+#' data(barents)
+#' Y <- barents$Abundance
+#' X <- model.matrix(Abundance ~ Latitude + Longitude + Depth + Temperature, data = barents)
+#' O <- log(barents$Offset)
+#' w <-- rep(1, nrow(Y))
+#' compute_PLN_starting_point(Y, X, O, w)
+#' }
+#'
+#' @importFrom stats lm.fit
+#' @export
+compute_PLN_starting_point <- function(Y, X, O, w, s = 0.1) {
+  # Y = responses, X = covariates, O = offsets (in log scale), w = weights
+  n <- nrow(Y); p <- ncol(Y); d <- ncol(X)
+  fits <- lm.fit(w * X, w * log((1 + Y)/exp(O)))
+  list(B = matrix(fits$coefficients, d, p),
+       M = matrix(fits$residuals, n, p),
+       S = matrix(s, n, p))
+}
diff --git a/_pkgdown.yml b/_pkgdown.yml
index 9cacd161..68f1a379 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -115,7 +115,7 @@ reference:
     - '`PLNfamily`'
     - '`plot.PLNfamily`'
     - '`rPLN`'
-    - '`rPLN`'
+    - '`compute_PLN_starting_point`'
 - title: Data sets
   desc: ~
   contents:
diff --git a/inst/check/jackknife_variance.R b/inst/check/jackknife_variance.R
new file mode 100644
index 00000000..a3fc2c05
--- /dev/null
+++ b/inst/check/jackknife_variance.R
@@ -0,0 +1,59 @@
+library(PLNmodels); library(mvtnorm)
+
+
+data("barents")
+n <- nrow(barents$Abundance); p <- ncol(barents$Abundance)
+fit <- PLN(Abundance ~ Latitude + Longitude + Depth + Temperature, data = barents)
+X <- model.matrix(Abundance ~ Latitude + Longitude + Depth + Temperature, data = barents)
+
+# Simulate data
+n <- 1e2
+X <- X[sample(1:nrow(X), n, replace=TRUE), ]
+B0 <- fit$model_par$B; Sigma0 <- fit$model_par$Sigma
+Y <- matrix(rpois(n*p, exp(X %*% B0 + mvtnorm::rmvnorm(n, sigma= Sigma0))), n, p)
+d <- ncol(X)
+
+# Jacknife sd estimate
+tictoc::tic()
+pln <- PLN(Y ~ -1 + X, control=PLN_param(config_post = list(jackknife=TRUE)))
+tictoc::toc()
+
+betaSdPackage <- as.vector(standard_error(pln, type='jackknife'))
+jackknife_estimates <- attr(pln$model_par$B, "estimates_jaccknife") |> purrr::map("B")
+betaJKPackage <- matrix(NA, n, p*d)
+for(i in 1:n){
+  betaJKPackage[i, ]  <- as.vector(jackknife_estimates[[i]])
+}
+
+
+# Jackknife sd estimate by hand
+tictoc::tic()
+betaJK <- matrix(NA, n, p*d)
+for(i in 1:n){
+  plnTmp <- PLN(Y[-i, ] ~ -1 + X[-i, ], control=PLN_param(trace=0))
+  betaJK[i, ]  <- as.vector(plnTmp$model_par$B)
+}
+tictoc::toc()
+betaSdByHand <- sqrt((n-1)*apply(betaJK, 2, function(x) { mean(x^2) - mean(x)^2 }))
+
+## Check that "automatic" and "manual" jackknife estimates are identical
+plot(as.vector(betaJK), as.vector(betaJKPackage),
+     ylab = "Computed from package", xlab = "Computed by hand",
+     main = "Jackknife estimates (all coefficients, all samples)")
+abline(0, 1, col = "red")
+max(abs(betaJK - betaJKPackage))
+
+# Empirical sd estimate
+B <- n
+betaEmp <- matrix(NA, B, p*d)
+for(b in 1:B){
+  if(b %% 10==0){cat(b, '')}
+  Ytmp <- matrix(rpois(n*p, exp(X %*% B0 + rmvnorm(n, sigma= Sigma0))), n, p)
+  plnTmp <- PLN(Ytmp ~ -1 + X, control=PLN_param(trace=0))
+  betaEmp[b, ]  <- as.vector(plnTmp$model_par$B)
+}
+cat("\n")
+betaSdEmpirical <- apply(betaEmp, 2, sd)
+
+# Comparison
+plot(as.data.frame(cbind(betaSdPackage, betaSdByHand, betaSdEmpirical)))
diff --git a/man/compute_PLN_starting_point.Rd b/man/compute_PLN_starting_point.Rd
new file mode 100644
index 00000000..03fb3ad7
--- /dev/null
+++ b/man/compute_PLN_starting_point.Rd
@@ -0,0 +1,39 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utils.R
+\name{compute_PLN_starting_point}
+\alias{compute_PLN_starting_point}
+\title{Helper function for PLN initialization.}
+\usage{
+compute_PLN_starting_point(Y, X, O, w, s = 0.1)
+}
+\arguments{
+\item{Y}{Response count matrix}
+
+\item{X}{Covariate matrix}
+
+\item{O}{Offset matrix (in log-scale)}
+
+\item{w}{Weight vector (defaults to 1)}
+
+\item{s}{Scale parameter for S (defaults to 0.1)}
+}
+\value{
+a named list of starting values for model parameter B and variational parameters M and S used in the iterative optimization algorithm of \code{\link[=PLN]{PLN()}}
+}
+\description{
+Barebone function to compute starting points for B, M and S when fitting a PLN. Mostly intended for internal use.
+}
+\details{
+The default strategy to estimate B and M is to fit a linear model with covariates \code{X} to the response count matrix (after adding a pseudocount of 1, scaling by the offset and taking the log). The regression matrix is used to initialize \code{B} and the residuals to initialize \code{M}. \code{S} is initialized as a constant conformable matrix with value \code{s}.
+}
+\examples{
+\dontrun{
+data(barents)
+Y <- barents$Abundance
+X <- model.matrix(Abundance ~ Latitude + Longitude + Depth + Temperature, data = barents)
+O <- log(barents$Offset)
+w <-- rep(1, nrow(Y))
+compute_PLN_starting_point(Y, X, O, w)
+}
+
+}
diff --git a/tests/testthat/test-standard-error.R b/tests/testthat/test-standard-error.R
index 9bd20803..6bdbf3aa 100644
--- a/tests/testthat/test-standard-error.R
+++ b/tests/testthat/test-standard-error.R
@@ -95,7 +95,7 @@ test_that("Check that variance estimation are coherent in PLNfit",  {
       trace           = 2
     )
 
-  myPLN$postTreatment(Y, X, exp(log_O), config = config_post)
+  myPLN$postTreatment(Y, X, log_O, config = config_post)
 
   tr_variational <- sum(standard_error(myPLN, "variational")^2)
   tr_bootstrap   <- sum(standard_error(myPLN, "bootstrap")^2)