results with centered and reduced W

vignettes/jSDM_Hmsc.Rmd

@@ -47,14 +47,21 @@ $$ \mathrm{g}(\theta_{ij}) = X_i\beta_j + W_i\lambda_j $$
 - $t_i$: number of visits at site $i$ 
 - $\mathrm{g}(\cdot)$: Link function (eg. logit or probit).
 
-The inference method is able to handle only one visit by site with a **probit** link function so $\forall i, \ t_i=1$ and $y_{ij} \sim \mathcal{B}ernoulli(\theta_{ij})$.
+The inference method is able to handle only one visit by site with a **probit** link function so $\forall i, \ t_i=1$ and $$y_{ij} \sim \mathcal{B}ernoulli(\theta_{ij})$$.  
 
-- $X_i$: Vector of explanatory variables for site $i$ (including intercept).
-- $\beta_j$: Effects of the explanatory variables on the probability of presence of species $j$.
-- $W_i$: Vector of random latent variables for site $i$. $W_i \sim N(0, 1)$. The number of latent variables must be fixed by the user (default to 2).
-- $\lambda_j$: Effects of the latent variables on the probability of presence of species $j$. Also known as "factor loadings" [@Warton2015].
+- $X_i$: Vector of explanatory variables for site $i$ with $X_i=(x_i^1,\ldots,x_i^p)\in \mathbb{R}^p$ where $p$ is the number of bio-climatic variables considered (including intercept $\forall i, x_i^1=1$).  
+
+- $\beta_j$: Effects of the explanatory variables on the probability of presence of species $j$ including species intercept ($\beta_{0j}$). We use a prior distribution $\beta_j \sim \mathcal{N}_p(0,\Sigma_{\beta})$ with $\Sigma_{\beta}$ a specified diagonal matrix of size $p \times p$.  
+
+- $W_i$: Vector of random latent variables for site $i$. $W_i \sim N(0, 1)$. The number of latent variables $q$ must be fixed by the user (default to $q=2$).  
+
+- $\lambda_j$: Effects of the latent variables on the probability of presence of species $j$ also known as "factor loadings" [@Warton2015]. We use the following prior distribution in the `jSDM` package to constraint values to $0$ on upper diagonal and to strictly positive values on diagonal, for $j=1,\ldots,J$ and $l=1,\ldots,q$ : $$\lambda_{jl} \sim \begin{cases}
+\mathcal{N}(0,1) & \text{if } l < j \\
+\mathcal{N}(0,1) \text{ left truncated by } 0 & \text{if } l=j \\
+P \text{ such as } \mathbb{P}(\lambda_{jl} = 0)=1  & \text{if } l>j
+\end{cases}$$.  
 
-This model is equivalent to a multivariate GLMM $\mathrm{g}(\theta_{ij}) = X_i.\beta_j + u_{ij}$, where $u_{ij} \sim \mathcal{N}(0, \Sigma)$ with the constraint that the variance-covariance matrix $\Sigma = \Lambda \Lambda^{\prime}$, where $\Lambda$ is the full matrix of factor loadings, with the $\lambda_j$ as its columns. 
+This model is equivalent to a multivariate GLMM $\mathrm{g}(\theta_{ij}) = X_i.\beta_j + u_{ij}$, where $u_{ij} \sim \mathcal{N}(0, \Sigma)$ with the constraint that the variance-covariance matrix $\Sigma = \Lambda \Lambda^{\prime}$, where $\Lambda$ is the full matrix of factor loadings, with the $\lambda_j$ as its columns.   
 
 Using this model we can compute the **full species residual correlation matrix** $R=(R_{ij})^{i=1,\ldots, nspecies}_{j=1,\ldots, nspecies}$ from the covariance in the latent variables such as : 
  $$\Sigma_{ij} = \begin{cases}
@@ -998,7 +1005,6 @@ save(T_jSDM_sim, Deviance_jSDM_sim, RMSE_jSDM_sim, RMSE_Hmsc_sim, T_Hmsc_sim, De
 ## Compilation time and deviance 
 
 ```{r time-deviance, echo=FALSE, eval=TRUE, out.width=900}
-
 load("jSDM_Hmsc_files/jSDM-Hmsc-comparison.rda")
 result <- data.frame(matrix(NA,5,7),row.names=c("Compilation time Hmsc (secondes)","Compilation time jSDM (secondes)", "Deviance Hmsc", "Deviance jSDM",""))
 colnames(result) <- c("Simulated","Mosquitos","Eucalypts","Frogs","Fungi","Birds","Mites") 
@@ -1033,7 +1039,6 @@ kable(result)%>%
 Computed for the probabilities of presences $\theta_{ij}$ on the simulated data-set.
 
 ```{r RMSE, eval=TRUE, echo=FALSE}
-
 result <- data.frame(matrix(NA,1,2),row.names=  c("RMSE"))
 colnames(result) <- c("Hmsc","jSDM")
 result$Hmsc <- RMSE_Hmsc_sim
@@ -1049,8 +1054,9 @@ kable(result, digits=3) %>%
 ### Simulated dataset
 
 ```{r jSDM-Hmsc-simulation, echo=FALSE}
+
 np <- ncol(X)
-print(paste(nrow(Y),"sites and ",ncol(Y)," species"),quote = FALSE)
+print(paste(nrow(Y),"sites and ", ncol(Y)," species"), quote = FALSE)
 # species fixed effect beta
 Hmsc_beta <- t(getPostEstimate(hM=mod_Hmsc_sim, parName='Beta')$mean)
 jSDM_beta <- matrix(0,nsp,np)

vignettes/jSDM_binomial_probit_sp_constrained.Rmd

@@ -50,36 +50,36 @@ library(kableExtra)
 
 Referring to the models used in the articles @Warton2015 and @Albert1993, we define the following model :
 
-$$ \mathrm{probit}(\theta_{ij}) =\alpha_i + \beta_{0j}+X_i.\beta_j+ W_i.\lambda_j $$
+$$ \mathrm{probit}(\theta_{ij}) =\alpha_i + X_i.\beta_j+ W_i.\lambda_j $$
+
+Considering the latent variable s$z$ such as $z_{ij} = \alpha_i + \beta_{0j} + X_i.\beta_j + W_i.\lambda_j + \epsilon_{i,j}$, with $\forall (i,j) \ \epsilon_{ij} \sim \mathcal{N}(0,1)$, it can be easily shown that: $y_{ij} \sim \mathcal{B}ernoulli(\theta_{ij})$, with: 
 
 - Link function probit: $\mathrm{probit}: q \rightarrow \Phi^{-1}(q)$ where $\Phi$ correspond to the distribution function of the reduced centered normal distribution.
 
-- Response variable: $Y=(y_{ij})^{i=1,\ldots,nsite}_{j=1,\ldots,nsp}$ with:
+- $\theta_{ij}$: occurrence probability of the species $j$ on site $i$ for $j=1,\ldots,J$ and for $i=1,\ldots,I$.
+
+- Response variable: $Y=(y_{ij})^{i=1,\ldots,nsite}_{j=1,\ldots,nsp}$ such as:
 
 $$y_{ij}=\begin{cases}
 0 & \text{ if species $j$ is absent on the site $i$}\\
 1 &  \text{ if species  $j$ is present on the site $i$}.
 \end{cases}$$
 
-- Latent variable $z_{ij} = \alpha_i + \beta_{0j} + X_i.\beta_j + W_i.\lambda_j + \epsilon_{i,j}$, with $\forall (i,j) \ \epsilon_{ij} \sim \mathcal{N}(0,1)$ and such that:
-
-$$y_{ij}=\begin{cases}
-1 & \text{if} \ z_{ij} > 0 \\
-0 &  \text{otherwise.}
-\end{cases}$$
-
-It can be easily shown that: $y_{ij} \sim \mathcal{B}ernoulli(\theta_{ij})$. 
+- $X_i$: Vector of explanatory variables for site $i$ with $X_i=(x_i^1,\ldots,x_i^p)\in \mathbb{R}^p$ where $p$ is the number of bio-climatic variables considered (including intercept $\forall i, x_i^1=1$).
 
-- Latent variables: $W_i=(W_i^1,\ldots,W_i^q)$ where $q$ is the number of latent variables considered, which has to be fixed by the user (by default $q=2$).
-We assume that $W_i \sim \mathcal{N}(0,I_q)$ and we define the associated coefficients: $\lambda_j=(\lambda_j^1,\ldots, \lambda_j^q)'$. We use a prior distribution $\mathcal{N}(0,10)$ for all lambdas not concerned by constraints to $0$ on upper diagonal and to strictly positive values on diagonal. 
+- $\alpha_i$: Random effect of site $i$ such as $\alpha_i \sim \mathcal{N}(0, V_{\alpha})$, corresponds to a mean suitability for site $i$. We assumed that $V_{\alpha} \sim \mathcal{IG}(\text{shape}=0.1, \text{rate}=0.1)$.
 
-- Explanatory variables: bioclimatic data about each site. $X=(X_i)_{i=1,\ldots,nsite}$ with $X_i=(x_i^1,\ldots,x_i^p)\in \mathbb{R}^p$ where $p$ is the number of bioclimatic variables considered.
-The corresponding regression coefficients for each species $j$ are noted : $\beta_j=(\beta_j^1,\ldots,\beta_j^p)'$.
+- $\beta_j$: Effects of the explanatory variables on the probability of presence of species $j$ including species intercept ($\beta_{0j}$). We use a prior distribution $\beta_j \sim \mathcal{N}_p(0,1)$. 
 
-- $\beta_{0j}$ correspond to the intercept for species $j$ which is assumed to be a fixed effect. We use a prior distribution $\mathcal{N}(0,10)$ for all betas. 
+- $W_i$: Vector of random latent variables for site $i$ such as $W_i \sim \mathcal{N}_q(0, 1)$ where the number of latent variables $q$ must be fixed by the user (default to $q=2$).
 
-- $\alpha_i$ represents the random effect of site $i$ such as $\alpha_i \sim \mathcal{N}(0,V_{\alpha})$ and we assume that $V_{\alpha} \sim \mathcal {IG}(\text{shape}=0.5, \text{rate}=0.005)$ as prior distribution by default. 
+- $\lambda_j$: Effects of the latent variables on the probability of presence of species $j$ also known as "factor loadings" [@Warton2015]. We use the following prior distribution to constraint values to $0$ on upper diagonal and to strictly positive values on diagonal, for $j=1,\ldots,J$ and $l=1,\ldots,q$ : $$\lambda_{jl} \sim \begin{cases}
+\mathcal{N}(0,1) & \text{if } l < j \\
+\mathcal{N}(0,1) \text{ left truncated by } 0 & \text{if } l=j \\
+P \text{ such as } \mathbb{P}(\lambda_{jl} = 0)=1  & \text{if } l>j
+\end{cases}$$.  
 
+This model is equivalent to a multivariate GLMM $\mathrm{g}(\theta_{ij}) = X_i.\beta_j + u_{ij}$, where $u_{ij} \sim \mathcal{N}(0, \Sigma)$ with the constraint that the variance-covariance matrix $\Sigma = \Lambda \Lambda^{\prime}$, where $\Lambda$ is the full matrix of factor loadings, with the $\lambda_j$ as its columns. 
 
 # Occurrence data-set
 
@@ -322,8 +322,8 @@ hist(mod_Fungi[[i]]$theta_latent,
 knitr::include_graphics(paste0("jSDM_binomial_probit_sp_constrained_files/figure-html/plot-results-probit-", 1:20, ".png"))
 ```
 
-Overall, the traces and the densities of the parameters indicate the convergence of the algorithm. Indeed, we observe on the traces that the values oscillate around averages without showing an upward or downward trend and we see that the densities are quite smooth and for the most part of Gaussian form, except for the factor loadings $lambda$.  
-In fact the values estimated for $lambda$ differ between the four chains and within each chain, the traces of these parameters oscillate very strongly showing upward or downward trends and we see that the densities are not of Gaussian form, as expected.  
+Overall, the traces and the densities of the parameters indicate the convergence of the algorithm. Indeed, we observe on the traces that the values oscillate around averages without showing an upward or downward trend and we see that the densities are quite smooth and for the most part of Gaussian form, except for the factor loadings $\lambda$.  
+In fact the values estimated for $\lambda$ differ between the four chains and within each chain, the traces of these parameters oscillate very strongly showing upward or downward trends and we see that the densities are not of Gaussian form, as expected.  
 
 This lack of convergence can be explained by the arbitrary choice of the species constrained to have positive values of $\lambda$ as being the first ones in the presence-absence data set. Indeed, the constrained species are supposed to be the ones that structure themselves most clearly on each latent axis but by arbitrarily choosing them this is not necessarily the case.
 
@@ -545,7 +545,7 @@ hist(mod_Fungi_ord[[i]]$theta_latent,
 knitr::include_graphics(paste0("jSDM_binomial_probit_sp_constrained_files/figure-html/plot-results-probit-ord-", 1:36, ".png"))
 ```
 
-Overall, the traces and the densities of the parameters indicate the convergence of the algorithm. Indeed, we observe on the traces that the values oscillate around averages without showing an upward or downward trend and we see that the densities are quite smooth and for the most part of Gaussian form, even for the factor loadings $lambda$ constrained to be positive.  
+Overall, the traces and the densities of the parameters indicate the convergence of the algorithm. Indeed, we observe on the traces that the values oscillate around averages without showing an upward or downward trend and we see that the densities are quite smooth and for the most part of Gaussian form, even for the factor loadings $\lambda$ constrained to be positive.  
 
 ### Matrice of residual correlations 
 

vignettes/jSDM_in_parallel.Rmd

@@ -70,7 +70,7 @@ library(jSDM)
 nsite <- 210
 
 #= Set seed for repeatability
-seed <- 123
+seed <- 1234
 set.seed(seed)
 
 #= Number of species
@@ -247,7 +247,7 @@ listData <- list(Y=Y, X=X[,-1], beta_start, lambda_start, W_start, alpha_start,
 
 
 ```{r mod-probit-2}
-
+ 
 ## Defining the function that will run MCMC on each CPU
 # Arguments:
 # i - will be 1 or 2
@@ -318,7 +318,8 @@ cat("content of each chain :", str_mod,"\n")
 
 The Gelman–Rubin diagnostic evaluates MCMC convergence by analyzing the difference between multiple Markov chains. The convergence is assessed by comparing the estimated between-chains and within-chain variances for each model parameter. Large differences between these variances indicate nonconvergence. See @Gelman1992 and @Brooks1998 for the detailed description of the method.
 
-Suppose we have $M$ chains, each of length $N$, although the chains may be of different lengths. The same-length assumption simplifies the formulas and is used for convenience. For a model parameter $\theta$, let $\left(\theta_{𝑚t}\right)_{t=1}^N$ be the $𝑚$th simulated chain, $𝑚=1,\dots,𝑀$. Let $\hat{\theta}_𝑚=\frac{1}{N}\sum\limits_{t=1}^N \hat{\theta}_{mt}$ and $\hat{\sigma}^2_𝑚=\frac{1}{N-1}\sum\limits_{t=1}^N (\hat{\theta}_{mt}-\hat{\theta}_𝑚)^2$ be the sample posterior mean and variance of the $𝑚$th chain, and let the overall sample posterior mean be $\hat{\theta}=\frac{1}{𝑀}\sum\limits_{m=1}^𝑀 \hat{\theta}_m$.
+Suppose we have $M$ chains, each of length $N$, although the chains may be of different lengths. The same-length assumption simplifies the formulas and is used for convenience. For a model parameter $\theta$, let $\left(\theta_{𝑚t}\right)_{t=1}^N$ be the $𝑚$th simulated chain, $𝑚=1,\dots,𝑀$.  
+Let $\hat{\theta}_𝑚=\frac{1}{N}\sum\limits_{t=1}^N \hat{\theta}_{mt}$ and $\hat{\sigma}^2_𝑚=\frac{1}{N-1}\sum\limits_{t=1}^N (\hat{\theta}_{mt}-\hat{\theta}_𝑚)^2$ be the sample posterior mean and variance of the $𝑚$th chain, and let the overall sample posterior mean be $\hat{\theta}=\frac{1}{𝑀}\sum\limits_{m=1}^𝑀 \hat{\theta}_m$.
 
 The between-chains and within-chain variances are given by
 $$B=\frac{N}{M-1}\sum\limits_{m=1}^𝑀 (\hat{\theta}_m - \hat{\theta})^2$$
@@ -389,7 +390,7 @@ Rhat <- data.frame(Rhat=c(psrf_alpha, psrf_V_alpha, psrf_beta, psrf_lambda, psrf
 library(ggplot2)
 ggplot(Rhat, aes(x=Variable, y=Rhat)) + 
   ggtitle("Averages of Rhat obtained for each type of parameter") +
-  theme(plot.title = element_text(hjust = 0.5, size=14)) +
+  theme(plot.title = element_text(hjust = 0.5, size=13)) +
   geom_bar(fill="skyblue", stat = "identity") +
   geom_text(aes(label=round(Rhat,3)), vjust=0, hjust=-0.1) +
   geom_hline(yintercept=1, color='red') +
@@ -402,38 +403,37 @@ knitr::include_graphics("jSDM_in_parallel_files/figure-html/MCMC-convergence-ran
 ```
 
 We can see that the $\hat{R}$ are very close to 1 for the species effects $\beta$, the site effects $\alpha$ and their variance $V_{alpha}$. We can therefore be fairly confident that convergence has been achieved for these parameters.    
+The $\hat{R}$ obtained for the latent variables $W$ and the factor loadings $\lambda$ are also very close to 1, indicating that convergence has been reached for these parameters.    
 
-However the $\hat{R}$ obtained for latent variables $W$ and factor loadings $\lambda$ are much greater than 1. It can be explained by the fact that $W$ and $\lambda$ can takes alternately positive and negative values between and even within MCMC chains, due to identifiability issues between $W$ and $\lambda$.  
-Consequently, if we consider the absolute values of the latent variables $W$ and  factor loadings $\lambda$ estimated to compute the $\hat{R}$, we get values well closer to 1. We can therefore consider that convergence has been almost reached for these parameters in absolute value. 
+However $\hat{R}$ associated to latent variables and factor loadings can be much greater than 1 in some cases. This can be explained by the fact that $W$ and $\lambda$ can alternatively take positive and negative values between and even within MCMC chains, due to identifiability issues between $W$ and $\lambda$.  
+Consequently, we also consider the absolute values of the latent variables $W$ and the factor loadings $\lambda$ estimated to compute the $\hat{R}$ in order to obtain values closer to 1. In this case, we can consider that convergence has been reached for these parameters in absolute value. 
 
 # Representation of results
 
 ```{r est-probit}
 ## Plot trace and posterior distributions 
 # for two first species
-plot(mcmc_list_param[,1:((np+n_latent)*2)])
+par(mfrow=c(2,2), cex.main=0.9)
+plot(mcmc_list_param[,1:((np+n_latent)*2)],
+     auto.layout=FALSE)
 # for two first sites 
-plot(mcmc_list_lv[,c(1:2,nsite+1:2)])
-par(mfrow=c(1,2))
+plot(mcmc_list_lv[,c(1:2,nsite+1:2)],
+     auto.layout=FALSE)
+par(mfrow=c(2,2))
 coda::traceplot(mcmc_list_V_alpha)
 coda::densplot(mcmc_list_V_alpha)
 abline(v=V_alpha.target, col='red')
 legend("topright", legend="V_alpha.target",
-       lwd=1,col='red', cex=0.6, bty="n")
-
-plot(mcmc_list_alpha[,c(1,2)])
+       lwd=1, col='red', cex=0.7, bty="n")
+plot(mcmc_list_alpha[,c(1,2)],
+     auto.layout=FALSE)
 # Deviance 
-plot(mcmc_list_deviance)
+plot(mcmc_list_deviance,
+     auto.layout=FALSE)
 ```
 
 ```{r est-probit-plot, echo=FALSE, out.width=800, eval=TRUE}
-knitr::include_graphics("jSDM_in_parallel_files/figure-html/est-probit-1.png")
-knitr::include_graphics("jSDM_in_parallel_files/figure-html/est-probit-2.png")
-knitr::include_graphics("jSDM_in_parallel_files/figure-html/est-probit-3.png")
-knitr::include_graphics("jSDM_in_parallel_files/figure-html/est-probit-4.png")
-knitr::include_graphics("jSDM_in_parallel_files/figure-html/est-probit-5.png")
-knitr::include_graphics("jSDM_in_parallel_files/figure-html/est-probit-6.png")
-knitr::include_graphics("jSDM_in_parallel_files/figure-html/est-probit-7.png")
+knitr::include_graphics(paste0("jSDM_in_parallel_files/figure-html/est-probit-", 1:9, ".png"))
 ```
 
 # Accuracy of predictions 
@@ -528,10 +528,7 @@ abline(a=0,b=1,col='red')
 ```
 
 ```{r obs-fitted-plot, echo=FALSE, out.width=800, eval=TRUE}
-knitr::include_graphics("jSDM_in_parallel_files/figure-html/obs-fitted-1.png")
-knitr::include_graphics("jSDM_in_parallel_files/figure-html/obs-fitted-2.png")
-knitr::include_graphics("jSDM_in_parallel_files/figure-html/obs-fitted-3.png")
-knitr::include_graphics("jSDM_in_parallel_files/figure-html/obs-fitted-4.png")
+knitr::include_graphics(paste0("jSDM_in_parallel_files/figure-html/obs-fitted-", 1:4, ".png"))
 ```
 
 # References