Hello,
This is a question in regards to Statistical Rethinking by McElreath(pg 517-522) cat example with missingness of a categorical variable (cat is present or not).
In the following stan model, the probability of cat is treated as a fixed effect:
mod_cat <- '
// stancode extracted from ulam object m15.31 (McElreath 2020, p. 519-520 )
data{
int N; // number of houses, where each house contains songbird
int RC[N]; // vector indicating unknown (missing value) whether cat in house (1 = missing, 0 if not)
int cat[N]; // vector indicating observed cat values (0,1,-9 if NA)
int notes[N]; // number of notes per house after sampling one minute
}
parameters{
real a; // intercept
real b; // parameter for association of cat presence/absence on number of notes
real<lower=0,upper=1> k; // parameter for probability that cat is present, not absent
}
model{
vector[N] lamda;
k ~ beta( 2 , 2); // prior for probability that cat is present
b ~ normal( 0 , 0.5 ); //prior for beta parameter
a ~ normal( 0 , 1 ); // prior for intercept
for ( i in 1:N ) {
lamda[i] = a + b * cat[i]; // rate of notes as function of cat present or absent
lamda[i] = exp(lamda[i]);
}
// likelihood of notes if cat not missing value (RC[i] = 0)
for ( i in 1:N )
if ( RC[i] == 0 ) cat[i] ~ bernoulli( k );
for ( i in 1:N )
if ( RC[i] == 0 ) notes[i] ~ poisson( lamda[i] );
// if cat is missing value (RC[i] = 1)
// average the probability of notes[i] for each possible value of cat (0,1)
// Pr(notes[i]) = Pr(cat[i] = 1)*(Pr(notes[i]|cat[i]=1) +
// Pr(cat[i]=0)*Pr(notes[i]|cat[i] = 0), where assume b = 0 because cat[i] = 0
// where Pr(cat[i] = 1) = k and Pr(cat[i] = 0) = (1-k)
for ( i in 1:N )
if ( RC[i] == 1 ) target += log_sum_exp(
log(k) + poisson_lpmf(notes[i] | exp(a + b)),
log(1 - k) + poisson_lpmf(notes[i] | exp(a))
);
}
}
'
I am trying to use this example for a different model but what if the cat being present or not was treated as a random effect as opposed to a fixed effect and here b is of size 2 instead of 1.
i.e.
mod_cat2 <- '
data{
int N; // number of houses, where each house contains songbird
int RC[N]; // vector indicating unknown (missing value) whether cat in house (1 = missing, 0 if not)
int cat[N]; // vector indicating observed cat values (0,1,-9 if NA)
int notes[N]; // number of notes per house after sampling one minute
int cat_miss_length; //here the categorical variable can only take values 0 or 1, so this should be set to '2'
}
parameters{
real a; // intercept
real b[cat_miss_length]; // parameter for association of cat presence/absence on number of notes
real<lower=0,upper=1> k; // parameter for probability that cat is present, not absent
}
model{
vector[N] lamda;
k ~ beta( 2 , 2); // prior for probability that cat is present
b ~ normal( 0 , 0.5 ); //prior for beta parameter
a ~ normal( 0 , 1 ); // prior for intercept
for ( i in 1:N ) {
lamda[i] = a + b * cat[i]; // rate of notes as function of cat present or absent
lamda[i] = exp(lamda[i]);
}
// likelihood of notes if cat not missing value (RC[i] = 0)
for ( i in 1:N )
if ( RC[i] == 0 ) cat[i] ~ bernoulli( k );
//following code is where there is a change 'b[1] and b[2]'
for ( i in 1:N )
if ( RC[i] == 0 ) notes[i] ~ poisson( lamda[i] );
for ( i in 1:N )
if ( RC[i] == 1 ) target += log_sum_exp(
log(k) + poisson_lpmf(notes[i] | exp(a + b[2])),
log(1 - k) + poisson_lpmf(notes[i] | exp(a+b[1]))
);
}
}
'
My understanding is this is fine to do since the likelihood specification works similarly in the fixed effect and random effect case. No source online says that this approach is fine to do and looking for confirmation to extend to a different model.
Thank you.