APPROACH

Variance based models and intraclass correlations (ICC) are approaches to examine the impact of so-called process variables on the measurements. This implementation is model-based.

Note: The term ICC is more frequently used to describe the agreement between different observers, examiners or even devices. In respective settings a good agreement is pursued. ICC-values can vary between \([-1; \: 1]\) and an ICC close to \(1\) is desired (Koo and Li 2016, Müller and Büttner 1994).

In multi-level analysis the ICC is interpreted differently. Please see Snijders et al. (Sniders and Bosker 1999). In this context the proportion of variance explained by respective group levels indicate an influence of (at least one) level of the respective group_vars.

Irrespective of the used terminology, regarding data quality it is desired that process variables do not explain systematically components of variance. Therefore, values close to \(0\) are desired.

ALGORITHM OF THIS IMPLEMENTATION:

  1. This implementation is yet restricted to data of type float.
  2. Missing codes are removed from resp_vars (if defined in the metadata)
  3. Deviations from limits, as defined in the metadata, are removed
  4. A linear mixed-effects model is estimated for resp_vars using co_vars and group_vars for adjustment.
  5. An output data frame is generated for group_vars indicating the ICC.

Example of study data

Data from the package dataquieR are loaded as shown below:

load(system.file("extdata", "study_data.RData", package = "dataquieR"))
sd1 <- study_data

This example of study data has N=3000 observations. Study data variables have abstract and non-interpretable names; appropriate labels must be mapped from the metadata.

Similarly to the approach of margins we assume that at least one examiner does not adhere to the SOP and may influence the measurement process.

v00000 v00001 v00002 v00003 v00004 v00005 v01003 v01002 v00103 v00006
3 LEIIX715 0 49 127 77 49 0 40-49 3.8
1 QHNKM456 0 47 114 76 47 0 40-49 1.9
1 HTAOB589 0 50 114 71 50 0 50-59 0.8
5 HNHFV585 0 48 120 65 48 0 40-49 3.8
1 UTDLS949 0 56 119 78 56 0 50-59 4.1
5 YQFGE692 1 47 133 81 47 1 40-49 9.5
1 AVAEH932 0 53 114 78 53 0 50-59 5.0
3 QDOPT378 1 48 116 86 48 1 40-49 9.6
3 BMOAK786 0 44 115 71 44 0 40-49 2.0
5 ZDKNF462 0 50 116 74 50 0 50-59 2.4

Example of metadata

Data from the package dataquieR are loaded as shown below:

load(system.file("extdata", "meta_data.RData", package = "dataquieR"))
md1 <- meta_data

Information corresponding to the study data is kept in the table of static metadata. An interpretable label for each variable is also attached.

Regarding the following implementation the columns DATA_TYPE, MISSING_LIST and HARD_LIMITS in the metadata are relevant. The R-function can be applied on variables of type float.

VAR_NAMES LABEL MISSING_LIST DATA_TYPE HARD_LIMITS
9 v00004 SBP_0 99980 | 99981 | 99982 | 99983 | 99984 | 99985 | 99986 | 99987 | 99988 | 99989 | 99990 | 99991 | 99992 | 99993 | 99994 | 99995 float [80;180]
10 v00005 DBP_0 99980 | 99981 | 99982 | 99983 | 99984 | 99985 | 99986 | 99987 | 99988 | 99989 | 99990 | 99991 | 99992 | 99993 | 99994 | 99995 float [50;Inf)
11 v00006 GLOBAL_HEALTH_VAS_0 99980 | 99983 | 99987 | 99988 | 99989 | 99990 | 99991 | 99992 | 99993 | 99994 | 99995 float [0;10]
14 v00009 ARM_CIRC_0 99980 | 99981 | 99982 | 99983 | 99984 | 99985 | 99986 | 99987 | 99988 | 99989 | 99990 | 99991 | 99992 | 99993 | 99994 | 99995 float [0;Inf)
21 v00014 CRP_0 99980 | 99981 | 99982 | 99983 | 99984 | 99985 | 99986 | 99988 | 99989 | 99990 | 99991 | 99992 | 99994 | 99995 float [0;Inf)
22 v00015 BSG_0 99980 | 99981 | 99982 | 99983 | 99984 | 99985 | 99986 | 99988 | 99989 | 99990 | 99991 | 99992 | 99994 | 99995 float [0;100]

Required R-packages

The above mentioned R-function requires several R-packages and utility functions:

library(ggplot2)

R-FUNCTION

The R-function requires the definition of up to nine arguments:

  • resp_vars: the name of one or more continuous measurement variables
  • group_vars: the name of one or more variables for the respective observer, device or reader. Caveat: the sequence of variables must correspond to the sequence of resp_vars
  • co_vars: a vector of covariables, e.g. age and sex for adjustment
  • label_col: the column of the metadata containing the label of variables
  • min_obs_in_subgroup: optional argument if a “group_var” is used. This argument specifies the minimum no. of observations that is required to include a subgroup (level) of the “group_var” in the analysis. Subgroups with less observations are excluded. The default is 30.
  • min_subgroups: optional argument if a “group_var” is used. This argument specifies the minimum no. of subgroups (levels) included “group_var.” If the variable defined in “group_var” has less subgroups it is not used for analysis. The default is 5.
  • study_data: the name of the data frame that contains the measurements
  • meta_data: the name of the data frame that contains metadata attributes of study data
acc_varcomp <-  function(resp_vars = NULL, group_vars, co_vars = NULL,
           min_obs_in_subgroup = 30, min_subgroups = 5, label_col = NULL,
           threshold_value = 0.05, study_data, meta_data) {

  # map meta to study
  util_prepare_dataframes()

  .min_obs_in_subgroup <- suppressWarnings(as.integer(min_obs_in_subgroup))
  if (is.na(.min_obs_in_subgroup)) {
   util_warning(c(
    "Could not convert min_obs_in_subgroup %s to a number.",
    "Set to standard value."
    ),
     dQuote(as.character(min_obs_in_subgroup)),
    applicability_problem = TRUE
   )
    min_obs_in_subgroup <- 30
  } else {
    min_obs_in_subgroup <- .min_obs_in_subgroup
  }

  .min_subgroups <- suppressWarnings(as.integer(min_subgroups))
  if (is.na(.min_subgroups)) {
    util_warning(
      "Could not convert min_subgroups %s to a number. Set to standard value.",
      dQuote(as.character(min_subgroups)),
      applicability_problem = TRUE
    )
    min_subgroups <- 5
  } else {
    min_subgroups <- .min_subgroups
  }

  # correct variable use?
  util_correct_variable_use("resp_vars",
    allow_more_than_one = TRUE,
    allow_null = TRUE,
    need_type = "integer|float"
  )

  util_correct_variable_use("group_vars",
    allow_null = TRUE,
    allow_any_obs_na = TRUE,
    allow_more_than_one = TRUE,
    need_type = "!float"
  )

  util_correct_variable_use("co_vars",
    allow_na = TRUE,
    allow_more_than_one = TRUE,
    allow_null = TRUE
  )

  rvs <- resp_vars

  if (length(rvs) == 0) {
    rvs <- colnames(ds1[, vapply(ds1, is.numeric, FUN.VALUE = logical(1))])
    rvs <- setdiff(rvs, group_vars)
    rvs <- setdiff(rvs, co_vars)
    if (length(group_vars) != 1)
      util_error(
        "Need exactly 1 group_vars, if all applicable resp_vars are used",
        applicability_problem = TRUE)
  }

  if (length(group_vars) == 1 && length(rvs) > 1) {
    group_vars <- rep(group_vars, length(rvs))
    message(sprintf("using the same group var %s for all resp_vars",
                    dQuote(group_vars[[1]])
                    ))
  }

  if (length(group_vars) != length(rvs)) {
    util_error(
      c("acc_varcomp expects one group_var per resp_var. Here,",
        "it has been called with %d resp_vars but %d group_vars"),
      length(rvs), length(group_vars),
      applicability_problem = TRUE)
  }

  ds1[, unique(group_vars)] <- lapply(ds1[, unique(group_vars), drop = FALSE],
                                      factor)

  # (3) if no covariables are defined for adjustment only the intercept is
  #     modelled
  if (is.null(co_vars) || length(co_vars) == 0) {
    co_vars <- "1"
  }

  # missing checks
  # -> encoded missings or jumps
  # -> resp_variables
  # -> variable types (factor(), ...)


  icc_output <- apply(cbind(rvs, group_vars), 1, function(row) {
    rv <- row[[1]]
    group_var <- row[[2]]

    check_df <- data.frame(table(ds1[[group_var]]))

    # too few observations in >1 level of the random effect
    critical_levels <- levels(check_df$Var1)[check_df$Freq <
                                               min_obs_in_subgroup]
    if (length(critical_levels) > 0) {
      util_warning("Levels %s were excluded due to less than %d observations.",
                   paste0(vapply(critical_levels, dQuote, ""), collapse = ", "),
                   min_obs_in_subgroup,
                   applicability_problem = FALSE)
      # exclude levels with too less observations
      ds1 <- ds1[!(ds1[[group_var]] %in% critical_levels), ]
      # dropping unused levels
      ds1[[group_var]] <- factor(ds1[[group_var]])
    }

    check_df <-
      data.frame(table(ds1[[group_var]])) # number of observers may have changed

    # less than min_subgroups levels of random effect
    if (length(check_df[, 1]) < min_subgroups) {
      util_warning("%d < %d levels in %s. Will not compute ICCs for %s.",
                   length(check_df[, 1]),
                   min_subgroups,
                   dQuote(group_var),
                   dQuote(rv),
                   applicability_problem = FALSE)
      return(data.frame(
        Variables = NA,
        Object = NA,
        Model.Call = NA,
        ICC = NA,
        Class.Number = NA,
        Mean.Class.Size = NA,
        Median.Class.Size = NA,
        Min.Class.Size = NA,
        Max.Class.Size = NA,
        convergence.problem = FALSE
      )[FALSE, , FALSE])
    }

    ds1[, rv] <- lapply(ds1[, rv, drop = FALSE], util_as_numeric)

    # prepare data frame for output
    res_df <- data.frame(
      Variables = NA,
      Object = NA,
      Model.Call = NA,
      ICC = NA,
      Class.Number = nrow(check_df),
      Mean.Class.Size = mean(check_df$Freq),
      Median.Class.Size = median(check_df$Freq),
      Min.Class.Size = min(check_df$Freq),
      Max.Class.Size = max(check_df$Freq),
      convergence.problem = FALSE
    )

    # fill output df
    res_df[1, 1] <- rv
    res_df[1, 2] <- group_var
    # create model formula from function arguments
    fmla <- as.formula(paste0(paste0(rv, "~"),
                              paste0(c(co_vars, paste0("(1|", group_var, ")")),
                                     collapse = "+")))
    # save formula for results
    res_df[1, 3] <- Reduce(paste, deparse(fmla))

    # remove misisng in used variables
    subdf <- na.omit(ds1[, c(rv, co_vars[co_vars %in% colnames(ds1)],
                             group_var), drop = FALSE])
    # drop unused levels (https://stackoverflow.com/a/1197154)
    subdf[] <- lapply(subdf, function(x) if (is.factor(x)) factor(x) else x)

    env <- environment()

    convergence.problem <- FALSE

    suppressMessages(withCallingHandlers(
      # fit mixed effects model with formula
      fit_lmer <- try(lme4::lmer(fmla, data = subdf, REML = TRUE)),
      message = function(m) {
        if (inherits(m, "simpleMessage") && any(
            grepl("singular",
             trimws(conditionMessage(m))))) {
          env$convergence.problem <- TRUE
        } else {
          message(m) # nocov
          # This is to handle other possible problems with the model
          # Since it is for unexpected problems, I cannot write a test
        }
      },
      warning = function(w) {
        # nocov start
        # on travis, the try above prevents a stop on convergence problems
        # for some unclear reasons. therefore, on travis,
        # the following code cannot be reached. However, locally, this
        # is possible.
        haystack <- trimws(paste0(conditionMessage(w), collapse = "\n"))
        if (inherits(w, "simpleWarning") && any(vapply(
            c("unable to evaluate scaled gradient",
              "Problem with Hessian check",
              "Model is nearly unidentifiable: very large eigenvalue",
              "convergence code [0-9\\-].+ from nloptwrap",
              "Model failed to converge with"
              ), grepl, haystack, perl = TRUE, FUN.VALUE = logical(1)))) {
          env$convergence.problem <- TRUE
          invokeRestart("muffleWarning")
        } else {
          warning(w)
          # This is to handle other possible problems with the model
          # Since it is for unexpected problems, I cannot write a test
        }
        # nocov end
      }
    ))

    res_df[1, "convergence.problem"] <- convergence.problem

    if (inherits(fit_lmer, "try-error")) {
      res_df[1, "convergence.problem"] <- TRUE
    } else {
      # extract variance of random effects
      v_tab <- as.data.frame(lme4::VarCorr(fit_lmer))
      res_df$Class.Number[[1]] <- lme4::ngrps(fit_lmer)
      # calculate icc
      res_df[1, 4] <- round(v_tab$vcov[1] / (v_tab$vcov[1] + v_tab$vcov[2]), 3)
    }

    class_sizes <- as.data.frame(table(subdf[, group_var]))
    res_df <- within(res_df, {
      Mean.Class.Size[[1]] <- mean(class_sizes$Freq)
      Median.Class.Size[[1]] <- median(class_sizes$Freq)
      Min.Class.Size[[1]] <- min(class_sizes$Freq)
      Max.Class.Size[[1]] <- max(class_sizes$Freq)
    })

    # output
    return(res_df)
  })

  icc_output <- as.data.frame(do.call(dplyr::bind_rows, icc_output),
                              stringsAsFactors = FALSE)
  rownames(icc_output) <- NULL

  max_icc <- suppressWarnings(max(icc_output$ICC, na.rm = TRUE))
  argmax_icc <- icc_output$Variables[which.max(icc_output$ICC)]

  icc_output[["GRADING"]] <- as.numeric(threshold_value <= icc_output$ICC)

  ## Format

  rownames(icc_output) <- NULL

  return(list(SummaryTable = icc_output, ScalarValue_max_icc = max_icc,
              ScalarValue_argmax_icc = argmax_icc))
}

Implementation and use of thresholds

This implementation makes no use of a threshold.

Call of the function

The following arguments have to be specified:

icc1 <- acc_varcomp(resp_vars = c("SBP_0", "DBP_0"), 
                group_vars = c("USR_BP_0", "USR_BP_0"),
                co_vars = c("AGE_0", "SEX_0"),
                label_col = "LABEL",
                min_obs_in_subgroup = 20,
                min_subgroups = 3,
                study_data = sd1,
                meta_data = md1)
names(icc1)
## [1] "SummaryTable"           "ScalarValue_max_icc"    "ScalarValue_argmax_icc"

OUTPUT

Output: Summary table

The summary data frame is called using icc1$SummaryTable:

Variables Object Model.Call ICC Class.Number Mean.Class.Size Median.Class.Size Min.Class.Size Max.Class.Size convergence.problem GRADING
SBP_0 USR_BP_0 SBP_0 ~ AGE_0 + SEX_0 + (1 | USR_BP_0) 0.153 15 165.8 160 29 413 FALSE 1
DBP_0 USR_BP_0 DBP_0 ~ AGE_0 + SEX_0 + (1 | USR_BP_0) 0.172 15 165.0 162 28 413 FALSE 1

In addition to this table some scalar values are returned (“ScalarValue_max_icc,” “ScalarValue_argmax_icc”) which represent the highest proportion ICC/VC and the response variable with the highest ICC/VC.

INTERPRETATION

ICC or the analysis of variance components should be applied in combination with MARGINS. Extended tests showed that ICC is less susceptible to false-positive indications of data quality issues than margins.

Limitations

Sufficient numbers of observations within each level of the group_vars are required. This can be specified by the formal min_obs_level. Nevertheless, the algorithm of the linear mixed effects model may not converge in cases of imbalanced and low numbers of observations.

Concept relations

Koo, T.K., and Li, M.Y. (2016). A guideline of selecting and reporting intraclass correlation coefficients for reliability research. Journal of Chiropractic Medicine 15, 155–163.
Müller, R., and Büttner, P. (1994). A critical discussion of intraclass correlation coefficients. Statistics in Medicine 13, 2465–2476.
Sniders, T., and Bosker, R. (1999). Multilevel analysis: An introduction to basic and advanced multilevel modeling. (Sage-Publications).