Calculate variable in generated quantities based on the maximum value of a vector

In the generated quantities block, the Stan program I am using has a 2-dimensional matrix Youden_index which is calculated for multiple diagnostic tests and thresholds.

For each test, I want to calculate the sensitivity and specificity at the threshold (the optimal threshold) which maximises the variable Youden_index. In other words, I need to create two vectors Se_at_opt and Sp_at_opt which correspond to the Se and Sp for each test calculated at their optimal thresholds.

Here is the code snippet in generated quantities:

generated quantities { 
                 for (t in 1:total_tests) { 
                     for (c in 1:num_categories[t]) { 
                        // Youden index
                        Youden_index[t,c] =  Se[t,c] + Sp[t,c] - 1;


            for (t in 1:total_tests) {  // t = test index
                 for (t_c in 1:num_thresholds[t]) { // t_c = threshold index

                     if  (Youden_index[t,t_c] ==  max(Youden_index[t, ])) {
                            Se_at_opt[t] = Se[t, t_c]; 
                            Sp_at_opt[t] = Sp[t, t_c]; 

However the outputs for Se_at_opt and Sp_at_opt are clearly wrong, since they do not match any of the outputs for the sensitivities and specificities calculated at all of the thresholds (i.e., Se_at_opt and Sp_at_opt do not mach any of the values for Se and Sp ).

It would be appreciated if somebody could advice as it seems like this should be a trivial thing to achieve, and I would rather have it done within the Stan program in generated quantities block rather than after the fact in R.

I can post full program and some dummy data if this would help


Did you manage to resolve this?

One possibility is that Youden_index[t,t_c] == max(Youden_index[t, ]) is always false, because there are some NaN values in Youden_index.

Best of luck with your model!

1 Like