greta_clique

greta_clique.R

# I am translating the code from https://github.com/mrcbarbier/diffuseclique/blob/master/codes/Clique%20Coexistence%20Matrix%20Generation.R to greta

# using the original function 
MatrixGen=function(etai, beta.m, beta.sd)
{
    S=length(etai)
    mat=diag(S)
    for(sp in 1:S)
    { 
        Nmi=etai[-sp]
        rn=matrix(rnorm((S-2)^2),nrow=S-2, ncol=S-2)
        tmp=rbind(Nmi[-1],rn)
        AA=rbind(Nmi,t(tmp))
        QQ=t(qr.Q(qr(AA)))
        QQ=QQ*sign(QQ[1,1])
        x=QQ%*%rep(beta.m,S-1)
        x[1]=(1-etai[sp])/sqrt(sum(Nmi^2))
        x[2:(S-1)]=x[2:(S-1)]+beta.sd*rnorm(S-2)
        row=t(QQ)%*%x
        mat[sp,-sp]=row
    }
    return(mat)
}

n_sp    <- 5
beta.m  <- 0.5
beta.sd <- 0.4

# the original authors chose arbitrary values of etai
# but I'm gonna simulate them from logits
(etai = plogis(rnorm(n_sp, -2, 1)))
m1 = MatrixGen(etai, beta.m, beta.sd)
image(m1)

# this could be how we generate N random interaction coef matrices
# I will call them "prior distributions"
m2 <- replicate(100, MatrixGen(etai, beta.m, beta.sd))

# this could be how we look at the correlations between priors
library(tidyverse)
m3 <- 
    reshape2::melt(m2) %>% 
    filter(Var1 != Var2) %>% 
    pivot_wider(names_from = c(Var1, Var2),
                values_from = value)
pairs(m3[, -1])


# now we code it in greta
library(greta)

S <- length(etai)
mat <- zeros(S, S)
diag(mat) <- 1
for (sp in 1:S) {
    Nmi = etai[-sp]
    rn = normal(0, 1, dim = c(S - 2, S - 2))
    tmp = rbind(t(Nmi[-1]), rn)
    AA = rbind(t(Nmi), t(tmp))
    
    # QQ = t(qr.Q(qr(AA)))
    r <- chol(t(AA) %*% AA)
    QQ <- t(AA %*% solve(r))
    
    QQ = QQ * sign(QQ[1, 1])
    x = QQ %*% greta_array(beta.m, c(S - 1, 1))
    x[1] = (1 - etai[sp]) / sqrt(sum(Nmi ^ 2))
    x[2:(S - 1)] = x[2:(S - 1)] + beta.sd * normal(0, 1, dim = S - 2)
    row = t(QQ) %*% x
    mat[sp, -sp] = row
}

# prior prediction using greta
mat_sim <- calculate(mat, nsim = 100)
dim(mat_sim$mat)
mat_sim_flat <- apply(mat_sim$mat, 1, as.vector)

# again looking at the prior correlations
library(bayesplot)
library(posterior)

mcmc_pairs(as_draws(t(mat_sim_flat)))

m4 <- 
    reshape2::melt(mat_sim$mat) %>% 
    filter(Var2 != Var3) %>% 
    pivot_wider(names_from = c(Var2, Var3),
                values_from = value)
pairs(m4[, -1])