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