`icc` and `get_variance` do not work inside custom functions #757

roaldarbol · 2024-07-22T13:59:47Z

roaldarbol
Jul 22, 2024

Okay, this is a super strange one! I've written a function that gather different ICC estimates... however, the performance and insight functions ( icc, get_variance) behave differently when inside the function than when they are called outside despite being specified exactly the same! I'm also using the rptR package, and have no issues with that. I can also use variance_decomposition with a Bayesian model without any issue. Any clue what might be going on here?

I found out that it fails if I fit the model inside the function, but it's works to specify the model outside the function and feed that as a parameter. I think it's super strange, and I have no idea why this is happening.

If you think it might be a problem with performance, then feel free to transfer it to an issue.

library(glmmTMB)
library(performance)
library(insight)
library(tibble)
library(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
set.seed(124)
d <- tibble::tibble(
  animal_id = factor(
    rep(
      c(
        "208", "209", "210", "214", "223", "228", "211", "213", "217", "222", "234",
        "241", "216", "230", "231", "240", "242", "244", "218", "220", "225", "237",
        "239", "219", "251", "252", "253", "254"
      ),
      each = 2L
    ),
    levels = c(
      "200", "204", "205", "206", "215", "224", "208", "209", "210", "214", "223",
      "228", "211", "213", "217", "222", "234", "241", "216", "230", "231", "240",
      "242", "244", "218", "220", "225", "237", "239", "245", "219", "236", "251",
      "252", "253", "254"
    )
  ),
  trial = rep(c(1, 2), 28),
  activity_ratio = c(
    0.09498024195328074, 0.06501857725091675, 0.11815590175658126,
    0.07432374123432427, 0.20160866416653947, 0.16168961981848032,
    0.18431356919353978, 0.1729636058911192, 0.08144358847160854,
    0.03429214493815157, 0.10657404730975248, 0.06867883079629465,
    0.05921577649099549, 0.10027697162542965, 0.22948745366059534,
    0.1318564983755886, 0.11521536711005413, 0.09988251038099903,
    0.11269870982289175, 0.10065372202090746, 0.30180913110275365,
    0.2784384586997261, 0.31090363146026273, 0.201433511508709,
    0.2064625890910772, 0.24142003722960942, 0.07379095189435501,
    0.05629768825555659, 0.23164588826907823, 0.2598283402933376,
    0.0773854788008177, 0.08506098926928175, 0.05543770346039902,
    0.08545932516134344, 0.13262226256298723, 0.13835676441921468,
    0.011223979573538618, 0.12175578019653631, 0.18844306077612688,
    0.18665264740941406, 0.30528203811103827, 0.36889354953270814,
    0.08101054157438016, 0.09939626397157338, 0.23026277596112824,
    0.3138068338604861, 0.10959072828854538, 0.1190276194587162,
    0.20927572336281877, 0.1531936170056966, 0.21758853563344488,
    0.12740603572231024, 0.11646060319509406, 0.05741400972805333,
    0.1767355830736662, 0.08745443933512403
  ),
) |>
  dplyr::group_by(animal_id)

beta_formula <- activity_ratio ~ trial + (1 | animal_id)
beta_model <- glmmTMB::glmmTMB(
  formula = as.formula(beta_formula),
  data = d,
  family = glmmTMB::beta_family()
)
insight::get_variance(beta_model)
#> $var.fixed
#> [1] 0.002488934
#> 
#> $var.random
#> [1] 0.3247858
#> 
#> $var.residual
#> [1] 0.08161196
#> 
#> $var.distribution
#> [1] 0.08161196
#> 
#> $var.dispersion
#> [1] 0
#> 
#> $var.intercept
#> animal_id 
#> 0.3247858
performance::icc(beta_model)
#> # Intraclass Correlation Coefficient
#> 
#>     Adjusted ICC: 0.799
#>   Unadjusted ICC: 0.794

# Here's the function we use to calculate all the ICCs for all tables
create_icc_summary <- function(
    df,
    beta_formula,
    ci_icc = NULL,
    ci_method = NULL
){
  print(list(ci_icc = ci_icc, ci_method = ci_method))
  
  # Get Beta ICC
  beta_model <- glmmTMB::glmmTMB(
    formula = as.formula(beta_formula),
    data = df,
    family = glmmTMB::beta_family()
  )
  
  insight::get_variance(beta_model)
  beta_icc <- performance::icc(beta_model)
  
  return(beta_icc)
}

create_icc_summary(
  df = d,
  beta_formula = activity_ratio ~ trial + (1 | animal_id),
  ci_icc = 0.95,
  ci_method = "analytical"
)
#> $ci_icc
#> [1] 0.95
#> 
#> $ci_method
#> [1] "analytical"
#> Warning: Can't calculate model's distribution-specific variance. Results are not
#>   reliable.
#>   A reason can be that the null model could not be computed manually. Try
#>   to fit the null model manually and pass it to `null_model`.
#> Warning: Can't calculate model's distribution-specific variance. Results are not
#>   reliable.
#>   A reason can be that the null model could not be computed manually. Try
#>   to fit the null model manually and pass it to `null_model`.
#> # Intraclass Correlation Coefficient
#> 
#>     Adjusted ICC: 1.000
#>   Unadjusted ICC: 0.992

# With the parameters pre-specified
create_icc_summary2 <- function(
    df,
    beta_formula
){
  beta_model <- glmmTMB::glmmTMB(
    formula = as.formula(beta_formula),
    data = df,
    family = glmmTMB::beta_family()
  )
  
  insight::get_variance(beta_model)
  performance::icc(beta_model)
}

create_icc_summary2(
  df = d,
  beta_formula = activity_ratio ~ trial + (1 | animal_id)
)
#> Warning: Can't calculate model's distribution-specific variance. Results are not
#>   reliable.
#>   A reason can be that the null model could not be computed manually. Try
#>   to fit the null model manually and pass it to `null_model`.
#> Warning: Can't calculate model's distribution-specific variance. Results are not
#>   reliable.
#>   A reason can be that the null model could not be computed manually. Try
#>   to fit the null model manually and pass it to `null_model`.
#> # Intraclass Correlation Coefficient
#> 
#>     Adjusted ICC: 1.000
#>   Unadjusted ICC: 0.992

# With the model given in advance
create_icc_summary_3 <- function(
    model
){

  insight::get_variance(model)
  performance::icc(model)
}

beta_model <- glmmTMB::glmmTMB(
  formula = as.formula(beta_formula),
  data = d,
  family = glmmTMB::beta_family()
)
create_icc_summary_3(
  beta_model
)
#> # Intraclass Correlation Coefficient
#> 
#>     Adjusted ICC: 0.799
#>   Unadjusted ICC: 0.794

^{Created on 2024-07-22 with reprex v2.1.0}

roaldarbol · 2024-07-22T14:02:34Z

roaldarbol
Jul 22, 2024
Author

Even if the model is specified outside, but it's specified inside a different function (but still available to the global environment), it doesn't work.

library(glmmTMB)
library(performance)
library(insight)
library(tibble)
library(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
set.seed(124)
d <- tibble::tibble(
  animal_id = factor(
    rep(
      c(
        "208", "209", "210", "214", "223", "228", "211", "213", "217", "222", "234",
        "241", "216", "230", "231", "240", "242", "244", "218", "220", "225", "237",
        "239", "219", "251", "252", "253", "254"
      ),
      each = 2L
    ),
    levels = c(
      "200", "204", "205", "206", "215", "224", "208", "209", "210", "214", "223",
      "228", "211", "213", "217", "222", "234", "241", "216", "230", "231", "240",
      "242", "244", "218", "220", "225", "237", "239", "245", "219", "236", "251",
      "252", "253", "254"
    )
  ),
  trial = rep(c(1, 2), 28),
  activity_ratio = c(
    0.09498024195328074, 0.06501857725091675, 0.11815590175658126,
    0.07432374123432427, 0.20160866416653947, 0.16168961981848032,
    0.18431356919353978, 0.1729636058911192, 0.08144358847160854,
    0.03429214493815157, 0.10657404730975248, 0.06867883079629465,
    0.05921577649099549, 0.10027697162542965, 0.22948745366059534,
    0.1318564983755886, 0.11521536711005413, 0.09988251038099903,
    0.11269870982289175, 0.10065372202090746, 0.30180913110275365,
    0.2784384586997261, 0.31090363146026273, 0.201433511508709,
    0.2064625890910772, 0.24142003722960942, 0.07379095189435501,
    0.05629768825555659, 0.23164588826907823, 0.2598283402933376,
    0.0773854788008177, 0.08506098926928175, 0.05543770346039902,
    0.08545932516134344, 0.13262226256298723, 0.13835676441921468,
    0.011223979573538618, 0.12175578019653631, 0.18844306077612688,
    0.18665264740941406, 0.30528203811103827, 0.36889354953270814,
    0.08101054157438016, 0.09939626397157338, 0.23026277596112824,
    0.3138068338604861, 0.10959072828854538, 0.1190276194587162,
    0.20927572336281877, 0.1531936170056966, 0.21758853563344488,
    0.12740603572231024, 0.11646060319509406, 0.05741400972805333,
    0.1767355830736662, 0.08745443933512403
  ),
) |>
  dplyr::group_by(animal_id)

# With the model given in advance
create_icc_summary_3 <- function(
    model
){

  insight::get_variance(model)
  performance::icc(model)
}

# Just a simple wrapper around glmmTMB
fit_model <- function(formula, df){
  glmmTMB::glmmTMB(
    formula = as.formula(formula),
    data = df,
    family = glmmTMB::beta_family()
  )
}

beta_formula <- activity_ratio ~ trial + (1 | animal_id)
mod <- fit_model(
  formula = beta_formula,
  df = d
)

create_icc_summary_3(
  model = mod
)
#> Warning: Can't calculate model's distribution-specific variance. Results are not
#>   reliable.
#>   A reason can be that the null model could not be computed manually. Try
#>   to fit the null model manually and pass it to `null_model`.
#> Warning: Can't calculate model's distribution-specific variance. Results are not
#>   reliable.
#>   A reason can be that the null model could not be computed manually. Try
#>   to fit the null model manually and pass it to `null_model`.
#> # Intraclass Correlation Coefficient
#> 
#>     Adjusted ICC: 1.000
#>   Unadjusted ICC: 0.992

beta_model <- glmmTMB::glmmTMB(
  formula = as.formula(beta_formula),
  data = d,
  family = glmmTMB::beta_family()
)
create_icc_summary_3(
  model = beta_model
)
#> # Intraclass Correlation Coefficient
#> 
#>     Adjusted ICC: 0.799
#>   Unadjusted ICC: 0.794

^{Created on 2024-07-22 with reprex v2.1.0}

0 replies

strengejacke · 2024-07-22T15:16:41Z

strengejacke
Jul 22, 2024
Maintainer

Could be a scoping issue. What if you fit the null-model and use the null_model argument?

1 reply

roaldarbol Jul 22, 2024
Author

Yeah, I imagine it's a scoping issue of some sort - I just have no clue what or how. null_model didn't do anything unfortunately.

library(glmmTMB)
library(performance)
library(insight)
library(tibble)
library(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
set.seed(124)
d <- tibble::tibble(
  animal_id = factor(
    rep(
      c(
        "208", "209", "210", "214", "223", "228", "211", "213", "217", "222", "234",
        "241", "216", "230", "231", "240", "242", "244", "218", "220", "225", "237",
        "239", "219", "251", "252", "253", "254"
      ),
      each = 2L
    ),
    levels = c(
      "200", "204", "205", "206", "215", "224", "208", "209", "210", "214", "223",
      "228", "211", "213", "217", "222", "234", "241", "216", "230", "231", "240",
      "242", "244", "218", "220", "225", "237", "239", "245", "219", "236", "251",
      "252", "253", "254"
    )
  ),
  trial = rep(c(1, 2), 28),
  activity_ratio = c(
    0.09498024195328074, 0.06501857725091675, 0.11815590175658126,
    0.07432374123432427, 0.20160866416653947, 0.16168961981848032,
    0.18431356919353978, 0.1729636058911192, 0.08144358847160854,
    0.03429214493815157, 0.10657404730975248, 0.06867883079629465,
    0.05921577649099549, 0.10027697162542965, 0.22948745366059534,
    0.1318564983755886, 0.11521536711005413, 0.09988251038099903,
    0.11269870982289175, 0.10065372202090746, 0.30180913110275365,
    0.2784384586997261, 0.31090363146026273, 0.201433511508709,
    0.2064625890910772, 0.24142003722960942, 0.07379095189435501,
    0.05629768825555659, 0.23164588826907823, 0.2598283402933376,
    0.0773854788008177, 0.08506098926928175, 0.05543770346039902,
    0.08545932516134344, 0.13262226256298723, 0.13835676441921468,
    0.011223979573538618, 0.12175578019653631, 0.18844306077612688,
    0.18665264740941406, 0.30528203811103827, 0.36889354953270814,
    0.08101054157438016, 0.09939626397157338, 0.23026277596112824,
    0.3138068338604861, 0.10959072828854538, 0.1190276194587162,
    0.20927572336281877, 0.1531936170056966, 0.21758853563344488,
    0.12740603572231024, 0.11646060319509406, 0.05741400972805333,
    0.1767355830736662, 0.08745443933512403
  ),
) |>
  dplyr::group_by(animal_id)

# With the model given in advance
create_icc_summary_3 <- function(
    model,
    null_mod
){
  
  insight::get_variance(model, null_model = null_mod)
  performance::icc(model, null_model = null_mod)
}

# Just a simple wrapper around glmmTMB
fit_model <- function(formula, df){
  glmmTMB::glmmTMB(
    formula = as.formula(formula),
    data = df,
    family = glmmTMB::beta_family()
  )
}

beta_formula <- activity_ratio ~ trial + (1 | animal_id)
mod <- fit_model(
  formula = beta_formula,
  df = d
)
null_mod <- null_model(mod)

create_icc_summary_3(
  model = mod,
  null_mod
)
#> Warning: Can't calculate model's distribution-specific variance. Results are not
#>   reliable.
#>   A reason can be that the null model could not be computed manually. Try
#>   to fit the null model manually and pass it to `null_model`.
#> Warning: Can't calculate model's distribution-specific variance. Results are not
#>   reliable.
#>   A reason can be that the null model could not be computed manually. Try
#>   to fit the null model manually and pass it to `null_model`.
#> # Intraclass Correlation Coefficient
#> 
#>     Adjusted ICC: 1.000
#>   Unadjusted ICC: 0.992

beta_model <- glmmTMB::glmmTMB(
  formula = as.formula(beta_formula),
  data = d,
  family = glmmTMB::beta_family()
)
create_icc_summary_3(
  model = beta_model
)
#> # Intraclass Correlation Coefficient
#> 
#>     Adjusted ICC: 0.799
#>   Unadjusted ICC: 0.794

^{Created on 2024-07-22 with reprex v2.1.0}

roaldarbol · 2024-07-23T08:15:05Z

roaldarbol
Jul 23, 2024
Author

I think there might be an issue that needs to be looked at - I've just knitted a document and again, the ICCs (especially the CIs) are off. When I run them line by line it works just fine. So whatever is happening it makes it un-reproducible. I'll upload a minimal project with data and .Rmd later today.

0 replies

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

Uh oh!

`icc` and `get_variance` do not work inside custom functions #757

Uh oh!

{{title}}

Uh oh!

Replies: 3 comments 1 reply

Uh oh!

{{title}}

Uh oh!

Uh oh!

{{title}}

Uh oh!

Uh oh!

{{title}}

Uh oh!

Uh oh!

{{title}}

Uh oh!

Uh oh!

{{editor}}'s edit

{{editor}}'s edit

Uh oh!

Select a reply

Uh oh!

Uh oh!

Uh oh!

icc and get_variance do not work inside custom functions #757

Uh oh!

roaldarbol Jul 22, 2024

Replies: 3 comments · 1 reply

Uh oh!

roaldarbol Jul 22, 2024 Author

Uh oh!

strengejacke Jul 22, 2024 Maintainer

Uh oh!

roaldarbol Jul 22, 2024 Author

Uh oh!

Uh oh!

roaldarbol Jul 23, 2024 Author

`icc` and `get_variance` do not work inside custom functions #757

roaldarbol
Jul 22, 2024

Replies: 3 comments 1 reply

roaldarbol
Jul 22, 2024
Author

strengejacke
Jul 22, 2024
Maintainer

roaldarbol Jul 22, 2024
Author

roaldarbol
Jul 23, 2024
Author