## Code used for the simulation study 
```{r second_sim_study, eval=F, echo=T, tidy=T, tidy.opts=list(width.cutoff=70)}
library(gllvm)
library(vegan)
R <- 1000
n <- 100
p <- 30
num.lv <- 2
set.seed(1)
beta0 <- matrix(runif(p,-1,1),ncol=p,nrow=n,byrow=T)
set.seed(1)
num.X <- 5
X <- mvtnorm::rmvnorm(n,rep(0,num.X),diag(num.X))
colnames(X)<-1:num.X
set.seed(1)
epsilon <- mvtnorm::rmvnorm(n,rep(0,num.lv),diag(num.lv))
epsilon <- resid(lm(epsilon~X)) #make sure that the residual is independent from the covariates as we assume
epsilon <- scale(epsilon, center = F)#get epsilon back to unit variance.
set.seed(1)
gamma<-matrix(runif(p*num.lv,-2,2),ncol=num.lv)
gamma[upper.tri(gamma)] <- 0
diag(gamma) <- 1
set.seed(1)
beta <- matrix(factanal(X,factors=2)$loadings,ncol=2)
qr.beta <- qr(beta)
beta <- qr.Q(qr.beta)%*%diag(diag(qr.R(qr.beta)))
LV <- (X%*%(beta)+epsilon)
eta <- beta0+LV%*%t(gamma)
eta3 <- beta0 + X%*%beta%*%t(gamma)
LVb <- X%*%beta

totresult <- NULL
  for(r in 1:R){
set.seed(r)
y <- matrix(rpois(p*n,exp(eta)),ncol=p,nrow=n)
set.seed(r)
y3 <- matrix(rbinom(p*n,size=1,prob=pnorm(eta)),ncol=p,nrow=n)
set.seed(r)
y5 <- matrix(rnorm(p*n,mean = eta),ncol=p,nrow=n)

#only fixed-effects
set.seed(r)
y7 <- matrix(rpois(p*n,exp(eta3)),ncol=p,nrow=n)
set.seed(r)
y8 <- matrix(rbinom(p*n,size=1,prob=pnorm(eta3)),ncol=p,nrow=n)
set.seed(r)
y9 <- matrix(rnorm(p*n,mean = eta3),ncol=p,nrow=n)

LV1 <- LV
LV3 <- LV
LV5 <- LV
LV7 <- LVb
LV8 <- LVb
LV9 <- LVb

#remove empty rows here for CA and CCA...
idx1 <- rowSums(y)
idx3 <- rowSums(y3)
idx7 <- rowSums(y7)
idx8 <- rowSums(y8)

X1 <- X[idx1>0,]
X3 <- X[idx3>0,]
X7 <- X[idx7>0,]
X8 <- X[idx8>0,]

y <- y[idx1>0,]
LV1 <- LV1[idx1>0,]
y3 <- y3[idx3>0,]
LV3 <- LV3[idx3>0,]
y7 <- y7[idx7>0,]
LV7 <- LV7[idx7>0,]
y8 <- y8[idx8>0,]
LV8 <- LV8[idx8>0,]

#Poisson
#Constrained simulation
modPC<-try(gllvm(y=y,X=X1,num.lv.c=num.lv,family="poisson",starting.val="res",maxit=1e7,num.lv=0, sd.errors=F, optimizer="alabama"),silent=T)
if(inherits(modPC,"try-error")){
modPC<-try(gllvm(y=y,X=X1,num.lv.c=num.lv,family="poisson",starting.val="res",maxit=1e7,num.lv=0, sd.errors=F, optimizer="alabama"),silent=T)
while(inherits(modPC,"try-error")){
modPC<-try(gllvm(y=y,X=X1,num.lv.c=num.lv,family="poisson",starting.val="res",maxit=1e7,num.lv=0, sd.errors=F, optimizer="alabama"),silent=T)
}
}
PCCA <- cca(y,X1)

#Fixed-effects simulation
modPCb<-try(gllvm(y=y7,X=X7,num.RR=num.lv,family="poisson",starting.val="zero",maxit=1e7,num.lv=0, sd.errors=F, optimizer="alabama"),silent=T)
if(inherits(modPCb,"try-error")){
modPCb<-try(gllvm(y=y7,X=X7,num.RR=num.lv,family="poisson",starting.val="res",maxit=1e7,num.lv=0, sd.errors=F, optimizer="alabama"),silent=T)
while(inherits(modPCb,"try-error")){
modPCb<-try(gllvm(y=y7,X=X7,num.RR=num.lv,family="poisson",starting.val="res",maxit=1e7,num.lv=0, sd.errors=F, optimizer="alabama"),silent=T)
}
}

PCCAb <- cca(y7,X7)

#Bernoulli
modBC<-try(gllvm(y=y3,X=X3,num.lv.c=num.lv,family="binomial",starting.val="res",maxit=1e7,num.lv=0, sd.errors=F, optimizer="alabama"),silent=T)

if(inherits(modBC,"try-error")){
modBC<-try(gllvm(y=y3,X=X3,num.lv.c=num.lv,family="binomial",starting.val="res",maxit=1e7,num.lv=0, sd.errors=F, optimizer="alabama"),silent=T)
while(inherits(modBC,"try-error")){
modBC<-try(gllvm(y=y3,X=X3,num.lv.c=num.lv,family="binomial",starting.val="res",maxit=1e7,num.lv=0, sd.errors=F, optimizer="alabama"),silent=T)
}
}else if(is.infinite(logLik(modBC))){
modBC<-try(gllvm(y=y3,X=X3,num.lv.c=num.lv,family="binomial",starting.val="res",maxit=1e7,num.lv=0, sd.errors=F, optimizer="alabama"),silent=T)
while(inherits(modBC,"try-error")){
modBC<-try(gllvm(y=y3,X=X3,num.lv.c=num.lv,family="binomial",starting.val="res",maxit=1e7,num.lv=0, sd.errors=F, optimizer="alabama"),silent=T)
}
while(is.infinite(logLik(modBC))){
modBC<-try(gllvm(y=y3,X=X3,num.lv.c=num.lv,family="binomial",starting.val="res",maxit=1e7,num.lv=0, sd.errors=F, optimizer="alabama"),silent=T)
}
}

BCCA <- cca(y3,X3)

#fixed-effects simulation
modBCb<-try(gllvm(y=y8,X=X8,num.RR=num.lv,family="binomial",starting.val="zero",maxit=1e7,num.lv=0, sd.errors=F, optimizer="alabama"),silent=T)
if(inherits(modBCb,"try-error")){
modBCb<-try(gllvm(y=y8,X=X8,num.RR=num.lv,family="binomial",starting.val="res",maxit=1e7,num.lv=0, sd.errors=F, optimizer="alabama"),silent=T)
while(inherits(modBCb,"try-error")){
modBCb<-try(gllvm(y=y8,X=X8,num.RR=num.lv,family="binomial",starting.val="res",maxit=1e7,num.lv=0, sd.errors=F, optimizer="alabama"),silent=T)
}
}else if(is.infinite(logLik(modBCb))){
modBCb<-try(gllvm(y=y8,X=X8,num.RR=num.lv,family="binomial",starting.val="res",maxit=1e7,num.lv=0, sd.errors=F, optimizer="alabama"),silent=T)
while(inherits(modBCb,"try-error")){
modBCb<-try(gllvm(y=y8,X=X8,num.RR=num.lv,family="binomial",starting.val="res",maxit=1e7,num.lv=0, sd.errors=F, optimizer="alabama"),silent=T)
}
while(is.infinite(logLik(modBCb))){
modBCb<-try(gllvm(y=y8,X=X8,num.RR=num.lv,family="binomial",starting.val="res",maxit=1e7,num.lv=0, sd.errors=F, optimizer="alabama"),silent=T)
}
}
BCCAb <- cca(y8,X8)

#Constrained simulation
modGC<-try(gllvm(y=y5,X=X,num.lv.c=num.lv,family="gaussian",starting.val="zero",maxit=1e7,num.lv=0, sd.errors=F, optimizer="alabama"),silent=T)
if(inherits(modGC,"try-error")){
modGC<-try(gllvm(y=y5,X=X,num.lv.c=num.lv,family="gaussian",starting.val="res",maxit=1e7,num.lv=0, sd.errors=F, optimizer="alabama"),silent=T)
while(inherits(modGC,"try-error")){
modGC<-try(gllvm(y=y5,X=X,num.lv.c=num.lv,family="gaussian",starting.val="res",maxit=1e7,num.lv=0, sd.errors=F, optimizer="alabama"),silent=T)
}
}else if(is.infinite(logLik(modGC))){
modGC<-try(gllvm(y=y5,X=X,num.lv.c=num.lv,family="gaussian",starting.val="res",maxit=1e7,num.lv=0, sd.errors=F, optimizer="alabama"),silent=T)
while(inherits(modGC,"try-error")){
modGC<-try(gllvm(y=y5,X=X,num.lv.c=num.lv,family="gaussian",starting.val="res",maxit=1e7,num.lv=0, sd.errors=F, optimizer="alabama"),silent=T)
}
while(is.infinite(logLik(modGC))){
modGC<-try(gllvm(y=y5,X=X,num.lv.c=num.lv,family="gaussian",starting.val="res",maxit=1e7,num.lv=0, sd.errors=F, optimizer="alabama"),silent=T)
}
}

RDA <- rda(y5,X)

#fixed-effects simulation
modGCb<-try(gllvm(y=y9,X=X,num.RR=num.lv,family="gaussian",starting.val="zero",maxit=1e7,num.lv=0, sd.errors=F, optimizer="alabama"),silent=T)
if(inherits(modGCb,"try-error")){
modGCb<-try(gllvm(y=y9,X=X,num.RR=num.lv,family="gaussian",starting.val="res",maxit=1e7,num.lv=0, sd.errors=F, optimizer="alabama"),silent=T)
while(inherits(modGCb,"try-error")){
modGCb<-try(gllvm(y=y9,X=X,num.RR=num.lv,family="gaussian",starting.val="res",maxit=1e7,num.lv=0, sd.errors=F, optimizer="alabama"),silent=T)
}
}else if(is.infinite(logLik(modGCb))){
modGCb<-try(gllvm(y=y9,X=X,num.RR=num.lv,family="gaussian",starting.val="res",maxit=1e7,num.lv=0, sd.errors=F, optimizer="alabama"),silent=T)
while(inherits(modGCb,"try-error")){
modGCb<-try(gllvm(y=y9,X=X,num.RR=num.lv,family="gaussian",starting.val="res",maxit=1e7,num.lv=0, sd.errors=F, optimizer="alabama"),silent=T)
}
while(is.infinite(logLik(modGCb))){
modGCb<-try(gllvm(y=y9,X=X,num.RR=num.lv,family="gaussian",starting.val="res",maxit=1e7,num.lv=0, sd.errors=F, optimizer="alabama"),silent=T)
}
}

RDAb <- rda(y9,X)

result <- matrix(0,ncol=5,nrow=12)
colnames(result) <- c("Simulation","Method","Distribution","SS","SS_spp")
result[,1] <- rep(c("Concurrent","Constrained"),times=6)
result[,2] <- c(rep(c("GLLVM","CCA"),each=4),rep(c("GLLVM","RDA"),each=2))
result[,3] <- c(rep(c(rep("Poisson",2),rep("Bernoulli",2)),times=2),rep("Gaussian",4))
result[,4] <- c(
procrustes(getLV(modPC),LV1,symmetric=T)$ss,
procrustes(getLV(modPCb),LV7,symmetric=T)$ss,
procrustes(getLV(modBC),LV3,symmetric=T)$ss,
procrustes(getLV(modBCb),LV8,symmetric=T)$ss,
procrustes(scores(PCCA,display="wa"),LV1,symmetric=T)$ss,
procrustes(scores(PCCAb, display="lc"),LV7,symmetric=T)$ss,
procrustes(scores(BCCA,display="wa"),LV3,symmetric=T)$ss,
procrustes(scores(BCCAb, display="lc"),LV8,symmetric=T)$ss,
procrustes(getLV(modGC,display="wa"),LV5,symmetric=T)$ss,
procrustes(getLV(modGCb, display="lc"),LV9,symmetric=T)$ss,
procrustes(scores(RDA,display="wa"),LV5,symmetric=T)$ss,
procrustes(scores(RDAb, display="lc"),LV9,symmetric=T)$ss
)
result[,5]<-c(
procrustes(modPC$params$theta,gamma,symmetric=T)$ss,
procrustes(modPCb$params$theta,gamma,symmetric=T)$ss,
procrustes(modBC$params$theta,gamma,symmetric=T)$ss,
procrustes(modBCb$params$theta,gamma,symmetric=T)$ss,
procrustes(scores(PCCA)$species,gamma,symmetric=T)$ss,
procrustes(scores(PCCAb)$species,gamma,symmetric=T)$ss,
procrustes(scores(BCCA)$species,gamma,symmetric=T)$ss,
procrustes(scores(BCCAb)$species,gamma,symmetric=T)$ss,
procrustes(modGC$params$theta,gamma,symmetric=T)$ss,
procrustes(modGCb$params$theta,gamma,symmetric=T)$ss,
procrustes(scores(RDA)$species,gamma,symmetric=T)$ss,
procrustes(scores(RDAb)$species,gamma,symmetric=T)$ss
)
totresult<-rbind(totresult,cbind(result,sim=r))
print(r)
  }
```