Estimate heritability from cumulative family model on an ordinal response?

I’m trying to estimate heritability of an ordinal response and I’m having a bit of trouble doing so. By big problem is that I don’t understand the output of cumulative family models and I don’ t understand where the residual variance is and how to put together the heritability formula from the brms model output. My data come from an experiment of plants grown in a common garden, and I have included in the model the row and column position of each genotype, as well as a random-intercept effect of genotype coded as “ind_code”.

For Gaussian family responses, I’ve been able to estimate heritability like so,

H2 = sd_ind_code__Intercept ^2 / (sd_ind_code__Intercept^2 + b_row^2 + b_col^2 + sigma^2)

However, for the ordinal data I haven’t been able to estimate reasonable values. Here’s the code I’ve been using for the

out <- posterior_samples(fit.dis1, pars = c("b_row", "b_col", "sd_ind_code__Intercept")) %>%
  mutate_at(.vars = 1:3, .funs = funs(.^2)) %>%
  mutate(H2 = sd_ind_code__Intercept / (sd_ind_code__Intercept + b_row + b_col )) %>%
  as_tibble()

data.frame(lwr5 = round(predictive_interval(as.matrix(out$H2))[1], 2),
           median = round(median(out$H2), 2),
           upr95 = round(predictive_interval(as.matrix(out$H2))[2], 2))
  lwr5 median upr95
     1      1     1

Now, It’s possible that heritability is 1, but I’ve estimated H2 in MCMCglmm and I know the value should be around 0.7.

The question I haven’t been able to answer is, where are the rest of the variance components from the model? Shouldn’t I be seeing a sigma in the model output that contains the residual variance?

Here’s a worked out example to recreate the problem with 100 samples.

dis1<- data.frame(
  ind_code=c("FIS_06","NEG_15","HWK_05","OFR_14","NEG_08","HBY_12","SLC_01","CPL_07","LPD_06","WTR_08","SSR_15","TUR_03","LPD_05","RMP_11","SSR_13","SLC_04","FIS_05","NEG_12","NEG_11","NEG_09","RMP_12","NEG_10","USDA10_08","SSR_12","DPR_07","NEG_06","DPR_06","USDA12_02","USDA12_03","CYH_13","TIM_12","TUR_04","LON_02","HBY_11","FIS_06","USDA13_04","OFR_11","USDA1_05","OUT_03","NBY_08","USDA4_07","USDA18_05","OUT_03","SLC_11","USDA12_04","NEG_03","WTR_12","OFR_05","JKH_04","CHL_06","RMP_09","NEG_08","HST_10","TIM_14","NEG_07","DCK_09","USDA15_04","WTR_10","USDA18_01","LON_02","MBK_14","CPL_11","FNO_14","LPD_04","USDA9_05","RAD_04","NEG_15","CLK_10","RMP_06","KAP_14","USDA19_04","SSR_02","USDA10_05","SSR_14","LPD_06","HBY_07","FIS_05","NEG_13","SKN_02","MMT_01","MBK_14","USDA12_10","CBI_05","USDA21_07","USDA5_10","OUT_10","USDA21_02","CYH_11","LSM_01","CHL_13","KAP_06","NBY_04","NEG_02","USDA4_06","MMT_11","USDA12_08","CLK_01","OFR_15","NEG_14","WTR_04"),
  rep_code=c("FIS_06_R1","NEG_15_R3","HWK_05_R4","OFR_14_R1","NEG_08_R1","HBY_12_R2","SLC_01_R2","CPL_07_R4","LPD_06_R3","WTR_08_R2","SSR_15_R2","TUR_03_R2","LPD_05_R3","RMP_11_R2","SSR_13_R2","SLC_04_R1","FIS_05_R1","NEG_12_R3","NEG_11_R3","NEG_09_R1","RMP_12_R3","NEG_10_R3","USDA10_08_R2","SSR_12_R2","DPR_07_R3","NEG_06_R2","DPR_06_R1.2","USDA12_02_R1","USDA12_03_R3","CYH_13_R3","TIM_12_R1","TUR_04_R4","LON_02_R2","HBY_11_R3","FIS_06_R2","USDA13_04_R1","OFR_11_R1","USDA1_05_R3","OUT_03_R4","NBY_08_R3","USDA4_07_R1","USDA18_05_R2","OUT_03_R5","SLC_11_R2","USDA12_04_R2","NEG_03_R1","WTR_12_R2","OFR_05_R3","JKH_04_R2","CHL_06_R2","RMP_09_R1","NEG_08_R2","HST_10_R2","TIM_14_R1","NEG_07_R1","DCK_09_R1.2","USDA15_04_R2","WTR_10_R1","USDA18_01_R1","LON_02_R1","MBK_14_R2","CPL_11_R4","FNO_14_R3","LPD_04_R3","USDA9_05_R2","RAD_04_R2","NEG_15_R1","CLK_10_R2","RMP_06_R3","KAP_14_R4","USDA19_04_R1","SSR_02_R3","USDA10_05_R3","SSR_14_R2","LPD_06_R2","HBY_07_R1","FIS_05_R2","NEG_13_R2","SKN_02_R4","MMT_01_R2","MBK_14_R4","USDA12_10_R1","CBI_05_R1","USDA21_07_R2","USDA5_10_R2","OUT_10_R4","USDA21_02_R1","CYH_11_R1","LSM_01_R1","CHL_13_R3","KAP_06_R5","NBY_04_R1","NEG_02_R2","USDA4_06_R3","MMT_11_R3","USDA12_08_R2","CLK_01_R2","OFR_15_R3","NEG_14_R2","WTR_04_R2"),
  d1=c(1,5,1,1,2,1,1,1,2,1,5,1,4,2,5,2,1,5,5,5,2,4,2,5,2,5,1,2,2,1,1,1,1,1,1,1,1,2,1,3,3,3,1,3,1,5,1,1,5,5,5,5,1,1,5,1,3,1,3,1,1,1,1,3,3,2,4,1,5,4,4,5,1,5,4,1,1,5,1,1,1,1,1,1,2,1,1,1,1,1,1,4,5,3,1,3,1,1,5,1),
  row=c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4),
  col=c(3,4,7,9,10,13,14,16,17,19,20,21,23,24,25,27,28,29,30,35,36,40,1,2,4,7,8,9,10,16,17,18,19,20,21,22,23,26,28,29,30,31,32,33,34,35,36,37,38,40,2,4,7,8,9,15,16,17,19,21,22,23,24,25,28,29,30,31,33,34,35,36,38,40,1,3,4,5,6,7,8,11,12,13,15,17,18,20,21,22,23,24,25,26,27,28,29,31,32,33))


fit.dis1 <- brm(d1 ~ row + col + (1|ind_code), 
                data = dis1,
                family = cumulative(threshold = "flexible"),
                chains = 2L, seed = 2398) %>% 
  add_criterion("loo", reloo = TRUE)


out <- posterior_samples(fit.dis1, pars = c("b_row", "b_col", "sd_ind_code__Intercept")) %>%
  mutate_at(.vars = 1:3, .funs = funs(.^2)) %>%
  mutate(H2 = sd_ind_code__Intercept / (sd_ind_code__Intercept + b_row + b_col )) %>%
  as_tibble()

data.frame(lwr5 = round(predictive_interval(as.matrix(out$H2))[1], 2),
           median = round(median(out$H2), 2),
           upr95 = round(predictive_interval(as.matrix(out$H2))[2], 2))

And here is what the model output of fit.dis1 looks like (using the full data set).

> fit.dis1
 Family: cumulative 
  Links: mu = logit; disc = identity 
Formula: d1 ~ row + col + (1 | ind_code) 
   Data: dis1 (Number of observations: 447) 
Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
         total post-warmup samples = 8000

Group-Level Effects: 
~ind_code (Number of levels: 282) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     3.25      0.46     2.44     4.24 1.00     1288     2451

Population-Level Effects: 
             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept[1]    -0.87      0.49    -1.88     0.09 1.00     3467     4776
Intercept[2]     0.39      0.49    -0.60     1.35 1.00     3436     4754
Intercept[3]     1.17      0.50     0.18     2.16 1.00     3256     4354
Intercept[4]     2.44      0.55     1.39     3.54 1.00     2807     3954
row             -0.07      0.02    -0.12    -0.02 1.00     3630     5468
col             -0.04      0.01    -0.07    -0.01 1.00     3768     5064

Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
>
> head(names(posterior_samples(fit.dis1)),10)
 [1] "b_Intercept[1]"               "b_Intercept[2]"              
 [3] "b_Intercept[3]"               "b_Intercept[4]"              
 [5] "b_row"                        "b_col"                       
 [7] "sd_ind_code__Intercept"       "r_ind_code[CBI_05,Intercept]"
 [9] "r_ind_code[CBI_08,Intercept]" "r_ind_code[CBI_10,Intercept]"

Can someone help me understand what I’m getting as output in an ordinal regression? What is the discrimination paramter in the model? Where is the residual variance (something akin to ‘sigma’ from gaussian regression models), and perhaps help me format the heritability calculation correctly?

Thanks all.

K

Version info:

  • Operating System: MacOS Catalina
  • R version: 4.0.0
  • brms Version: 2.12.0

I don’t know what heritability and such is, but if you’re having trouble figuring out what is going on inside brms you can use the brms functions make_stancode and make_standata to look at the actual data and model being passed in to Stan.

It it super helpful for situations like this to make sure you’re computing what you think you are.