-
Notifications
You must be signed in to change notification settings - Fork 1
/
measure_variances.R
95 lines (89 loc) · 4.04 KB
/
measure_variances.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
#' Measure variances from a covariance matrix
#'
#' Creates a table of variance components derived from a user-defined covariance matrix.
#'
#' @param mat A symmetric \code{n x n} positive (semi)-definite variance matrix.
#' @param estimate When \code{TRUE} (default is \code{FALSE}), the variance components are treated
#' as estimates and sample corrections are applied.
#'
#' @return A table which partitions the total variance into main effect and interaction variance,
#' heterogeneity of variance and lack of correlation, and non-crossover and crossover interaction.
#' When \code{correction = TRUE}, a sample correction is applied to the variance components.
#'
#' @examples
#' # Generate a structured covariance matrix and then measure the partitioning of variance.
#'
#' diag_mat <- rand_diag_mat(
#' n = 10
#' )
#' cor_mat <- struc_cor_mat(
#' n = 10,
#' )
#'
#' cov_mat <- sqrt(diag_mat) %*% cor_mat %*% sqrt(diag_mat)
#'
#' measure_variances(
#' mat = cov_mat
#' )
#'
#' @export
measure_variances <- function(mat,
estimate = FALSE) {
if (!is.matrix(mat)) stop("'mat' must be a matrix")
if (!isSymmetric(mat)) stop("'mat' must be a symmetric matrix")
if (any(diag(mat) < 0)) stop("All diagonal elements of 'mat' must >= 0")
if (any(is.na(mat))) stop("'mat' must not contain missing values")
if (!is.logical(estimate)) stop("'estimate' must be logical")
n <- ncol(mat)
total_var <- mean(diag(mat))
if (estimate) {
main_eff <- mean(mat[upper.tri(mat)])
if (main_eff < 0) stop("The mean upper triangular element of 'mat' must be >= 0")
if (main_eff < 1e-8) warning("The main effect variance is zero")
int_eff <- mean(diag(mat - main_eff))
if (round(main_eff + int_eff, 8) != round(total_var, 8)) stop("Sum of main effect and interaction variances does not match total variance")
sqrt_dg <- sqrt(diag(mat))
het_scale <- sum((sqrt_dg - mean(sqrt_dg))^2) / (n - 1)
lack_cor <- sum(outer(sqrt_dg, sqrt_dg, "*") * (1 - stats::cov2cor(mat))) / n / (n - 1)
if (round(het_scale + lack_cor, 8) != round(int_eff, 8)) stop("Sum of heterogeneity of scale and lack of correlation does not match interaction variance")
df <- data.frame(
Component = c(
"Main effect", "Interaction",
"Heterogeniety of scale", "Lack of correlation",
"Total"
),
Value = c(main_eff, int_eff, het_scale, lack_cor, total_var)
)
} else if (!estimate) {
main_eff <- mean(mat)
if (main_eff < 0) stop("The mean element of 'mat' must be >= 0")
if (main_eff < 1e-8) warning("The main effect variance is zero")
int_eff <- mean(diag(mat - main_eff))
if (round(main_eff + int_eff, 8) != round(total_var, 8)) stop("Sum of main effect and interaction variances does not match total variance")
sqrt_dg <- sqrt(diag(mat))
het_scale <- sum((sqrt_dg - mean(sqrt_dg))^2) / n
lack_cor <- sum(outer(sqrt_dg, sqrt_dg, "*") * (1 - stats::cov2cor(mat))) / n^2
if (round(het_scale + lack_cor, 8) != round(int_eff, 8)) stop("Sum of heterogeneity of scale and lack of correlation does not match interaction variance")
cov_mat <- rowMeans(mat)
adjust <- min(0, cov_mat)
if (adjust < 0) {
warning("Non-crossover and crossover interaction not completely independent, covariances will be attributed to crossover variance")
}
main_eff_adjust <- main_eff - adjust
cov_mat_adjust <- cov_mat - adjust
non_cross <- main_eff * t(cov_mat_adjust) %*% cov_mat_adjust / n / (main_eff_adjust^2)
cross <- total_var - main_eff * t(cov_mat_adjust) %*% cov_mat_adjust / n / (main_eff_adjust^2)
if (round(non_cross + cross, 8) != round(total_var, 8)) stop("Sum of non-crossover and crossover variance does not match total variance")
df <- data.frame(
Component = c(
"Main effect", "Interaction",
"Heterogeniety of scale", "Lack of correlation",
"Non-crossover", "Crossover",
"Total"
),
Value = c(main_eff, int_eff, het_scale, lack_cor, non_cross, cross, total_var)
)
}
df$Proportion <- df$Value / total_var
return(df)
}