Tag capture-recapture von Bertalanffy hierarchical model

Modeling
,
1

Hi guys, I would like to know how to estimate parameters Linf, t0 and k in an hierarchical model of growth von Bertalanffy http://www.pisces-conservation.com/growthhelp/index.html?von_bertalanffy.htm, incorporating the discrete variables “J” and “S”, given:

N: total of individuals
J: fishes that have been recaptured 1,2, …6 times
L: length at recapture (L)
dt: time between tagging and recapture
S: places (1 and 2)

cat(file = “VBN.stan”, "
data {
int<lower=1> N;
int S;
vector[N] L;
int<lower=1,upper=S> J[N];
}

parameters {               
real<lower=0> k;        
real<lower=40> Linf;  
real t0;    
real<lower=0> sd_k;
real<lower=0> sd_linf;
real<lower=0> sd_t0;
real<lower=0> sigma;      
real z_k[S];
real z_linf[S];
real z_t0[S];
vector[N] A;

}
transformed parameters {
vector[S] r_linf = sd_linf * z_linf;
vector[S] r_k = sd_k * z_k;
vector[S] r_t0 = sd_t0 * z_t0;
vector[N] mu = r_linf[S] * (1 - exp( - r_k[S] * (A - r_t0[S])));
}
model {
Linf ~ normal(50, 10);
k ~ normal(0.3, 0.15);
t0 ~ cauchy(0, 0.01);
sigma ~ cauchy(0, 5);
sd_linf ~ cauchy(0, 5);
z_linf ~ normal(0, 1);
sd_k ~ cauchy(0, 1);
z_k ~ cauchy(0, 1);
sd_t0 ~ cauchy(0, 0.01);
z_t0 ~ normal(0, 1);
A ~ normal(0, 0.1);
L ~ normal(mu, sigma);
}
generated quantities {
vector[N] Y_pred;
for (i in 1:N) {
Y_pred[i] = normal_rng(mu[i], sigma);
}
}
")
I have problems with dimensions declarated
I appreciate any advice, thanks!

2

The model will at least compile if you replace:

real z_k[S];
real z_linf[S];
real z_t0[S];

with:

vector[S] z_k;
vector[S] z_linf;
vector[S] z_t0;

Vectors (the second thing) and arrays (the first thing) behave differently in Stan. Only vectors support multiplication with scalars, which the model is using here:

vector[S] r_linf = sd_linf * z_linf;
vector[S] r_k = sd_k * z_k;
vector[S] r_t0 = sd_t0 * z_t0;

That help?

3

The model will at least compile if you replace:

real z_k[S];
real z_linf[S];
real z_t0[S];

with:

vector[S] z_k;
vector[S] z_linf;
vector[S] z_t0;

Vectors (the second thing) and arrays (the first thing) behave differently in Stan. Only vectors support multiplication with scalars, which the model is using here:

vector[S] r_linf = sd_linf * z_linf;
vector[S] r_k = sd_k * z_k;
vector[S] r_t0 = sd_t0 * z_t0;

That help?

I try this

cat(file = “VBN.stan”, "
data {
int N;
int J;
int S;
vector[N] L;
vector<lower = 0>[N] dt;
int<lower=1,upper=J> j[N];
int<lower=1,upper=S> s[N];
}

parameters {
vector<lower=0>[J] k;        
vector<lower=40>[S] Linf;  
real t0;    
real tau;      
vector<lower = 0>[N] A;   
real centre;
real precage;
real<lower=0> sd_k;
real<lower=0> sd_linf;
real z_k;
real z_linf;

}
transformed parameters {
vector[N] mu = Linf[s] .* (1 - exp( - k[j] .* (A - t0 + dt)));
real logcentre;
real sigma;
logcentre=log(centre);
sigma=1/(tau*tau);
}
model {
Linf ~ normal(sd_linf, z_linf);
k ~ normal(sd_k, z_k);
tau ~ gamma(.0001, .0001);
sigma ~ normal(0, 5);
centre ~ uniform(0,50);
precage ~ gamma(.01,.01);
A ~ lognormal(logcentre,precage);
t0 ~ cauchy(0, 0.01);
sd_linf ~ normal(50, 10);
z_linf ~ normal(0, 1);
sd_k ~ cauchy(0, 1);
z_k ~ cauchy(0, 1);
L ~ normal(mu, sigma);
}
generated quantities {
vector[N] Y_pred;
for (i in 1:N) {
Y_pred[i] = normal_rng(mu[i], sigma);
}
}
")

But, I have problems with dimensions declared , gave me

mismatch in number dimensions declared and found in context; processing stage=initialization; variable name=k; dims declared=(6); dims found=()"

My data is:
dato<-read.csv(“nmarca.csv”)
L<-dato$A1
dt<-dato$dt
N<-NROW(dato)
J<-length(unique(dato$marca))
S<-length(unique(dato$site))
j<-dato$marca
s<-dato$site

and in stan_data:

stan_data <- list(N, L, dt,S,J,s,j)
inits <- function() list(Linf = runif(1, 45, 55), k = runif(1, 0, 1),sigma = runif(1, 0, 1))
fit0 <- stan(file = “VBN.stan”, data = stan_data, init = inits, iter = 1000,
thin = 1, warmup=500, chains = 3)

I can’t fit effect to discrete variable (S and J) in the model, please, Anyone has any advice, I will be very thanks.

4

k is declared to be length J:

vector<lower=0>[J] k;

But you’re only passing in something of length 1 to the inits:

k = runif(1, 0, 1)

Assuming J is a variable declared to be the right thing, you probably want:

k = runif(J, 0, 1)

That work?

5

Yes, work, thanks.

I would like knows how fit effect of “J” in deviation “sigma”, I change the effect of “S” in:

vector[N] mu = Linf[s] .* (1 - exp( - k[s] .* (A - t0 + dt)));

But, How to do effect of “J” (numbers times that was measured the fish, 6 times) in desviation.

6

Do you mean estimating J itself? So estimating if J == 1 or if J == 2 or something?

Or do you mean estimating J different values of sigma? Like, if sigma were declared like:

real sigma[J];