Iterating in my model (ordered_logistic) gives wrong answer only on the first iteration

I am building a model that has data from subpopulations. I am gradually working up to partial pooling but at this point I just want to have separate estimates for each subpopulation and another for the combined subpopulations. I am making the model as explicit as possible. Optimization may or may not come later.

Below is a printout with the toy data I am using and the model and the results. Notice that the first subpopulation sp[1,] has estimated means of all .25’s (the prior means) far from the data, yet the second subpop has an estimates that match the data and the population estimates also seem correct. I suspect there is some kind of subscripting error. Any hints how I can check the stan code for this type of error before running it would also be appreciated. Thanks for any help.
n = n # total number of measurements calculated with sum(sum.counts)

print(n)
[1] 37
k = nsubpop # number of subpopulations calculated by lrn - frn + 1
print(k)
[1] 2
KK = samp.size # vector of subpopulation sample sizes
print(KK)
[1] 21 16
max_sample_size = max.samp.size # fixed maximum sample
print(max_sample_size)
[1] 40
all_rubric_scores = rubric_score # all population measurements combined
print(all_rubric_scores)
[1] 1 1 1 1 1 1 1 1 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4
sub_rubric_scores = sub_rubric_score # subpopulation measurements
print(sub_rubric_scores)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20]
[1,] 1 1 1 2 2 2 3 3 3 3 3 3 3 3 3 3 4 4 4 4
[2,] 1 1 1 1 1 2 2 2 3 3 3 4 4 4 4 4 100 100 100 100
[,21] [,22] [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37] [,38]
[1,] 4 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100
[2,] 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100
[,39] [,40]
[1,] 100 100
[2,] 100 100

set phi to 0. It is a placeholder for later correlations

phi = 0
tmu = .25 # prior mu
tsigma = .3 # prior sigma

stancode = "

  • data {
    
  •   int<lower=0> n; // total number of measurements
    
  •   int<lower=0> k; // number of subpopulations 
    
  •   int KK[k]; // vector of subpopulation sample sizes
    
  •   int<lower=0> max_sample_size; // fixed maximum subpopulation sample size
    
  •   int all_rubric_scores[n]; // all population measurements combined
    
  •   int sub_rubric_scores[k,max_sample_size] ; // subpopulation measurements
    
  •   real tmu;
    
  •   real tsigma;
    
  •   real phi; // placeholder
    
  •     }
    
  • parameters {
    
  •   simplex[4] p;  // probabilities for the total population. they add to 1
    
  •   simplex[4] sp[k]; // probabilities for subpopulations
    
  • }
    
  • transformed parameters {
    
  •   vector[3] cum_p;  // total population cumulative probabilities
    
  •   vector[3] cutpoints; // total population cutpoints
    
  •   matrix[3,k] cum_sp; // subpopulation cumulative probabilities
    
  •   matrix[3,k] sp_cutpoints; // subpop cutpoints
    
  •   cum_p = cumulative_sum(p[1:3]); // total pop - only three needed since cum_p[4] is always 1
    
  •   cutpoints = logit(cum_p);  // total pop -logit is the inverse logistic
    
  •   for (i in 1:k) {cum_sp[,i] = cumulative_sum(sp[i,1:3]); // subpop cumulative probabilities
    
  •                   sp_cutpoints[,i] = logit(cum_sp[,i]); // subpop cutpoints
    
  •                   };  
    
  •  }
    
  • model {
    
  •   // Total population
    
  •   for (i in 1:n) {
    
  •       all_rubric_scores[i]~ordered_logistic(phi,cutpoints);
    
  •                   }
    
  •   for (i in 1:4)  {
    
  •     p[i] ~normal(tmu, tsigma) T[0, 1]; // probability priors truncated normal
    
  •                    }
    
  •   // Sub populations
    
  •   for (i in 1:k)  {
    
  •           for(j in 1:KK[k])
    
  •                     { sub_rubric_scores[k,j]~ordered_logistic(phi,sp_cutpoints[,k]);
    
  •                     }
    
  •                   }
    
  •   for (i in 1:k)  {
    
  •           for (j in 1:4)  {
    
  •                             sp[i,j]~normal(tmu, tsigma) T[0, 1]; // probability priors truncated normal
    
  •                           }
    
  •                  }
    
  • }
    
  • "

writeLines(stancode,“mrubric.stan”)
data_list = list(n=n,k=k,KK=KK,all_rubric_scores=all_rubric_scores,max_sample_size=max_sample_size,sub_rubric_scores=sub_rubric_scores,phi=phi,tmu=tmu,tsigma=tsigma)
mrubricord = stan(“mrubric.stan”,data = data_list,chains=4,cores=4,iter = 2000)
Warning in system(paste(CXX, ARGS), ignore.stdout = TRUE, ignore.stderr = TRUE) :
‘-E’ not found

show(mrubricord)
Inference for Stan model: mrubric.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.

                 mean se_mean   sd    2.5%     25%     50%     75%   97.5% n_eff Rhat

p[1] 0.22 0.00 0.06 0.11 0.18 0.22 0.26 0.36 7446 1
p[2] 0.17 0.00 0.06 0.08 0.13 0.17 0.21 0.30 7835 1
p[3] 0.34 0.00 0.07 0.21 0.29 0.33 0.38 0.48 8127 1
p[4] 0.27 0.00 0.07 0.15 0.22 0.26 0.31 0.40 8356 1
sp[1,1] 0.25 0.00 0.16 0.01 0.12 0.23 0.36 0.61 5567 1
sp[1,2] 0.25 0.00 0.16 0.01 0.12 0.23 0.36 0.61 6132 1
sp[1,3] 0.25 0.00 0.16 0.01 0.12 0.23 0.36 0.61 6795 1
sp[1,4] 0.25 0.00 0.16 0.01 0.12 0.23 0.35 0.60 6374 1
sp[2,1] 0.30 0.00 0.07 0.17 0.25 0.30 0.35 0.45 7407 1
sp[2,2] 0.20 0.00 0.07 0.09 0.15 0.19 0.24 0.34 7367 1
sp[2,3] 0.20 0.00 0.06 0.09 0.15 0.19 0.24 0.34 6752 1
sp[2,4] 0.30 0.00 0.07 0.17 0.25 0.30 0.35 0.46 8160 1
cum_p[1] 0.22 0.00 0.06 0.11 0.18 0.22 0.26 0.36 7446 1
cum_p[2] 0.40 0.00 0.08 0.25 0.34 0.40 0.45 0.55 7709 1
cum_p[3] 0.73 0.00 0.07 0.60 0.69 0.74 0.78 0.85 8356 1
cutpoints[1] -1.29 0.00 0.38 -2.10 -1.54 -1.27 -1.03 -0.57 7192 1
cutpoints[2] -0.43 0.00 0.32 -1.08 -0.64 -0.43 -0.21 0.19 7566 1
cutpoints[3] 1.04 0.00 0.34 0.40 0.80 1.03 1.26 1.75 8322 1
cum_sp[1,1] 0.25 0.00 0.16 0.01 0.12 0.23 0.36 0.61 5567 1
cum_sp[1,2] 0.30 0.00 0.07 0.17 0.25 0.30 0.35 0.45 7407 1
cum_sp[2,1] 0.50 0.00 0.19 0.14 0.36 0.50 0.63 0.87 6051 1
cum_sp[2,2] 0.50 0.00 0.08 0.35 0.45 0.50 0.55 0.65 7369 1
cum_sp[3,1] 0.75 0.00 0.16 0.40 0.65 0.77 0.88 0.99 6374 1
cum_sp[3,2] 0.70 0.00 0.07 0.54 0.65 0.70 0.75 0.83 8160 1
sp_cutpoints[1,1] -1.41 0.02 1.24 -4.56 -2.00 -1.21 -0.60 0.46 3070 1
sp_cutpoints[1,2] -0.86 0.00 0.35 -1.59 -1.09 -0.85 -0.62 -0.21 6741 1
sp_cutpoints[2,1] 0.00 0.01 0.91 -1.83 -0.56 -0.01 0.55 1.89 5620 1
sp_cutpoints[2,2] 0.00 0.00 0.32 -0.62 -0.22 0.00 0.21 0.60 7344 1
sp_cutpoints[3,1] 1.38 0.02 1.16 -0.42 0.61 1.19 1.97 4.20 4519 1
sp_cutpoints[3,2] 0.87 0.00 0.36 0.18 0.62 0.86 1.10 1.59 8357 1
lp__ -112.45 0.06 2.35 -118.04 -113.79 -112.10 -110.74 -108.89 1597 1

Samples were drawn using NUTS(diag_e) at Thu Dec 16 13:56:12 2021.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).

1 Like

DONT BOTHER. It was a subscripting error!

1 Like