Broad-Sense Heritability using H2cal function from inti package

General
,
1

How can I generalize this function below using either purrr or for loop to calculate the broad sense heritability for more than one triats. The code will work for a single traits.
I have tried several step but I am not getting it.

install.packages("inti")

dt <- potato
hr <- H2cal(data = dt
, trait = "stemdw"
, gen.name = "geno"
, rep.n = 5
, fixed.model = "0 + (1|bloque) + geno"
, random.model = "1 + (1|bloque) + (1|geno)"
, emmeans = TRUE
, plot_diag = TRUE
, outliers.rm = TRUE
)

hr$model %>% summary()

hr$tabsmr %>% kable(caption = "Variance component table")

hr$blues %>% kable(caption = "BLUEs")

hr$blups %>% kable(caption = "BLUPs")

Thank you very much!!!

2

Here's how to create a list of results of running the model on each trait. Extracting components may be a bit of a pain.

library(inti)
#> Loading required package: shiny
#> Loading required package: ggplot2
#> Loading required package: dplyr
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
#> Loading required package: tidyr
#> Loading required package: tibble
library(gt)

dt <- potato
run_model <- function(x) {
  H2cal(data = dt
              , trait = x
              , gen.name = "geno"
              , rep.n = 5
              , fixed.model = "0 + (1|bloque) + geno"
              , random.model = "1 + (1|bloque) + (1|geno)"
              , emmeans = TRUE
              , plot_diag = TRUE
              , outliers.rm = TRUE
  )
}


traits <- colnames(potato)[4:17]
result <- sapply(traits,run_model)
#> boundary (singular) fit: see help('isSingular')
#> boundary (singular) fit: see help('isSingular')
#> boundary (singular) fit: see help('isSingular')
#> boundary (singular) fit: see help('isSingular')
#> boundary (singular) fit: see help('isSingular')


#> boundary (singular) fit: see help('isSingular')
#> boundary (singular) fit: see help('isSingular')
#> boundary (singular) fit: see help('isSingular')
#> boundary (singular) fit: see help('isSingular')


#> boundary (singular) fit: see help('isSingular')
#> boundary (singular) fit: see help('isSingular')
#> boundary (singular) fit: see help('isSingular')
#> boundary (singular) fit: see help('isSingular')


#> boundary (singular) fit: see help('isSingular')
#> boundary (singular) fit: see help('isSingular')
#> boundary (singular) fit: see help('isSingular')
#> boundary (singular) fit: see help('isSingular')


#> boundary (singular) fit: see help('isSingular')
#> boundary (singular) fit: see help('isSingular')
#> boundary (singular) fit: see help('isSingular')
#> boundary (singular) fit: see help('isSingular')


#> boundary (singular) fit: see help('isSingular')
#> boundary (singular) fit: see help('isSingular')
#> boundary (singular) fit: see help('isSingular')
#> boundary (singular) fit: see help('isSingular')



result[1]
#> [[1]]
#>     trait rep geno env year     mean      std   min   max      V.g      V.e
#> 1 spad_29   5   15   1    1 56.36487 4.137189 46.93 61.96 16.19783 9.185086
#>        V.p repeatability      H2.s      H2.p      H2.c
#> 1 18.03485     0.8981407 0.8981407 0.9463373 0.9463373
result[2]
#> [[1]]
#> # A tibble: 15 × 2
#>    geno  spad_29
#>    <chr>   <dbl>
#>  1 G01      55.1
#>  2 G02      53.2
#>  3 G03      57.2
#>  4 G04      47.4
#>  5 G05      55.9
#>  6 G06      60.3
#>  7 G07      58.7
#>  8 G08      61.7
#>  9 G09      51.2
#> 10 G10      55.5
#> 11 G11      61.1
#> 12 G12      54.3
#> 13 G13      59.0
#> 14 G14      59.7
#> 15 G15      55.1
result[3]
#> [[1]]
#> # A tibble: 15 × 6
#>    geno  spad_29    SE    df lower.CL upper.CL
#>    <fct>   <dbl> <dbl> <dbl>    <dbl>    <dbl>
#>  1 G01      55.1 0.958  129.     53.2     57.0
#>  2 G02      53.0 0.958  129.     51.1     54.9
#>  3 G03      57.2 0.958  129.     55.3     59.1
#>  4 G04      46.9 0.958  129.     45.0     48.8
#>  5 G05      55.9 0.958  129.     54.0     57.8
#>  6 G06      60.5 0.958  129.     58.6     62.4
#>  7 G07      58.8 0.958  129.     56.9     60.7
#>  8 G08      62.0 0.958  129.     60.1     63.9
#>  9 G09      50.9 0.958  129.     49.1     52.8
#> 10 G10      55.4 0.958  129.     53.5     57.3
#> 11 G11      61.4 0.958  129.     59.5     63.3
#> 12 G12      54.2 0.958  129.     52.3     56.1
#> 13 G13      59.2 0.958  129.     57.3     61.1
#> 14 G14      59.9 0.958  129.     58.0     61.8
#> 15 G15      55.0 0.958  129.     53.1     56.9
result[4]
#> [[1]]
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: spad_29 ~ 1 + (1 | bloque) + (1 | geno)
#>    Data: dt.rm
#> Weights: weights
#> REML criterion at convergence: 799.2244
#> Random effects:
#>  Groups   Name        Std.Dev.
#>  geno     (Intercept) 4.025   
#>  bloque   (Intercept) 0.000   
#>  Residual             3.031   
#> Number of obs: 150, groups:  geno, 15; bloque, 5
#> Fixed Effects:
#> (Intercept)  
#>       56.36  
#> optimizer (nloptwrap) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings
result[5]
#> [[1]]
#> [[1]]$fixed
#>  [1] bloque      geno        spad_29     resi        res_MAD     rawp.BHStud
#>  [7] index       adjp        bholm       out_flag   
#> <0 rows> (or 0-length row.names)
#> 
#> [[1]]$random
#>  [1] bloque      geno        spad_29     resi        res_MAD     rawp.BHStud
#>  [7] index       adjp        bholm       out_flag   
#> <0 rows> (or 0-length row.names)

result[7]
#> [[1]]
#>     trait rep geno env year     mean      std   min   max     V.g      V.e
#> 1 spad_83   5   15   1    1 42.01864 3.662827 37.14 50.59 11.5687 17.60003
#>       V.p repeatability      H2.s      H2.p      H2.c
#> 1 15.0887     0.7667126 0.7667126 0.8556439 0.8661048

Created on 2023-08-08 with reprex v2.0.2

3

Thank you very much for the help @technocrat, I finally stick to the purrr way which I used the code below:

Extension of the models for all the trait

library(tidyverse)
library(inti)
library(knitr)

dt <- potato

mod <- function(dat,name){
hr <- H2cal(data = dat
, trait = "value"
, gen.name = "geno"
, rep.n = 5
, fixed.model = "0 + (1|bloque) + geno"
, random.model = "1 + (1|bloque) + (1|geno)"
, emmeans = TRUE
, plot_diag = FALSE
, outliers.rm = TRUE
)
return(hr)
#hr$tabsmr %>% kable(caption = "Variance component table")
}

Nested data frame

dat <- dt |>
pivot_longer(cols = c(4:17)) |>
group_by(name) |>
nest() %>%
mutate(model=map2(.x = data,.y=name,.f = mod )) |>
select(-data)

dat %>%
pluck('model') %>%
set_names(dat$name)

4

This topic was automatically closed 21 days after the last reply. New replies are no longer allowed.

If you have a query related to it or one of the replies, start a new topic and refer back with a link.