Hello everyone!

I’m trying to practice phylogenetic regressions and I’d like to fit them in raw Stan, rather than brms or another package.

I’m working with R simulation code from Ives’ book on correlated data, but I’ve written the Stan model myself. I’m looking for feedback on the way I wrote this model. Should I be using a Cholesky decomposition somehow? is there a more flexible or general approach?

## simulation code

```
## simulate data
set.seed(1618)
n <- 20
b0 <- 0
b1 <- .5
lam.x <- .98
lam.e <- .5
phy <- ape::compute.brlen(
ape::rtree(n=n),
method = "Grafen",
power = 1)
plot(phy)
```

```
phy.x <- phylolm::transf.branch.lengths(
phy=phy, model="lambda",
parameters=list(lambda = lam.x))$tree
phy.e <- phylolm::transf.branch.lengths(
phy=phy, model="lambda",
parameters=list(lambda = lam.e))$tree
x <- ape::rTraitCont(phy.x, model = "BM", sigma = 1)
e <- ape::rTraitCont(phy.e, model = "BM", sigma = .3)
x <- x[match(names(e), names(x))]
Y <- b0 + b1 * x + e
Y <- array(Y)
rownames(Y) <- phy$tip.label
plot(x, Y)
```

^{Created on 2024-09-27 with reprex v2.1.1}

## Stan code:

```
data {
int n;
int s;
vector[n] x;
vector[n] y;
matrix[s, s] phyvcv;
}
parameters {
real b0;
real b1;
real sigma_x;
real sigma_y;
real logit_lambda_x;
real logit_lambda_y;
}
transformed parameters {
real<lower=0,upper=1> lambda_x;
lambda_x = inv_logit(logit_lambda_x);
// y
real<lower=0,upper=1> lambda_y;
lambda_y = inv_logit(logit_lambda_y);
}
model {
b0 ~ std_normal();
b1 ~ normal(.5, .5);
sigma_x ~ exponential(1);
sigma_y ~ exponential(1);
logit_lambda_x ~ normal(3, .2);
logit_lambda_y ~ normal(0, .2);
matrix[s, s] vcv_x;
vcv_x = add_diag(sigma_x^2*lambda_x*phyvcv, sigma_x^2*(1 - lambda_x));
matrix[s, s] vcv_y;
vcv_y = add_diag(sigma_y^2*lambda_y*phyvcv, sigma_y^2*(1 - lambda_y));
x ~ multi_normal(rep_vector(0, n), vcv_x);
y ~ multi_normal(b0 + b1*x, vcv_y);
}
```

## Results

Parameter recovery is pretty good! very encouraging, but once again I would love any constructive feedback! Thanks :)

```
> phylo_sample
variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
lp__ 23.73 24.08 1.90 1.71 20.07 26.12 1.00 1771 2248
b0 -0.08 -0.08 0.12 0.12 -0.28 0.12 1.00 2671 2639
b1 0.55 0.55 0.09 0.08 0.41 0.69 1.00 2482 2482
sigma_x 0.97 0.94 0.16 0.15 0.74 1.27 1.00 3658 2083
sigma_y 0.22 0.21 0.04 0.04 0.16 0.29 1.00 2989 1939
logit_lambda_x 3.01 3.01 0.20 0.19 2.68 3.33 1.00 4703 3031
logit_lambda_y 0.00 0.00 0.20 0.20 -0.33 0.34 1.00 4438 2770
lambda_x 0.95 0.95 0.01 0.01 0.94 0.97 1.00 4703 3031
lambda_y 0.50 0.50 0.05 0.05 0.42 0.58 1.00 4438 2770
```