Background

When building variance-covariance matrices from parameter estimation results, one problem might be that the resulting matrix is not positive definite. One solution to that problem might be a technique called bending. There are different approaches how bending can be implemented.

Different Approaches

An first intuitive approach is to decrease all off-diagonal elements by a small amount and to increase the diagonal elements also by some small quantity. This is done iteratively until the smallest eigenvalue of the resulting variance-covariance matrix is larger than some specified lower limit.

Schaeffer Method

Alternatively, the eigenvalues below a certain threshold can be projected above that threshold leaving the ratios of the distances between the eigenvalues constant. This is implemented in the function makePD2(). Hence for a non-positive definite matrix containing estimation results, the following steps can be used to bend the matrix.

(mat_npd <- matrix(data = c(100, 80, 20, 6, 80, 50, 10, 2, 20, 10, 6, 1, 6, 2, 1, 1), nrow = 4))
#>      [,1] [,2] [,3] [,4]
#> [1,]  100   80   20    6
#> [2,]   80   50   10    2
#> [3,]   20   10    6    1
#> [4,]    6    2    1    1

Based on the eigenvalues mat_npd is non-positive definite

eigen(mat_npd, only.values = TRUE)$values
#> [1] 162.1627196   4.1339019   0.9171925 -10.2138140

The matrix mat_npd can be bent using

(mat_bent1 <- makePD2(A = mat_npd))
#>            [,1]      [,2]      [,3]     [,4]
#> [1,] 103.615081 75.499246 18.377481 5.013141
#> [2,]  75.499246 55.603411 12.020026 3.228634
#> [3,]  18.377481 12.020026  6.728218 1.442922
#> [4,]   5.013141  3.228634  1.442922 1.269397

As can be seen, the eigenvalues of the bent matrix are all positive, but the ratio between the smallest and the largest eigenvalue is big.

eigen(mat_bent1, only.values = TRUE)$values
#> [1] 1.621627e+02 4.133902e+00 9.171925e-01 2.292926e-03

Ratio of Eigenvalues

If this ratio is of any concern, the second bending function make_pd_rat_ev() can be used which allows for the specification of a maximum ratio between smallest and largest eigenvalue.

(mat_bent2 <- make_pd_rat_ev(A = mat_npd, pn_max_ratio = 100))
#>            [,1]      [,2]       [,3]      [,4]
#> [1,] 104.195930 74.774864 18.0837820 4.9954987
#> [2,]  74.774864 56.506986 12.3914112 3.2288164
#> [3,]  18.083782 12.391411  7.0139582 0.8658603
#> [4,]   4.995499  3.228816  0.8658603 3.7720338

Leading to a much smaller range of eigenvalues as can be seen from the output below.

eigen(mat_bent2, only.values = TRUE)$values
#> [1] 162.162720   4.133902   3.570658   1.621627

Jorjani Method

In Jorjani et al. 2003, the authors described weighted and unweighted bending. The unweighted version allows the user to specify a minimum eigenvalue. All eigenvalues of the input matrix that are smaller are just set to this minimum. Then the result matrix is reconstructed using the eigenvector-eigenvalue decomposition with the modified eigenvalues. This can be done as follows

(mat_input <- matrix(data = c(100,95,80,40,40,
                                   95,100,95,80,40,
                                   80,95,100,95,80,
                                   40,80,95,100,95,
                                   40,40,80,95,100), ncol = 5))
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]  100   95   80   40   40
#> [2,]   95  100   95   80   40
#> [3,]   80   95  100   95   80
#> [4,]   40   80   95  100   95
#> [5,]   40   40   80   95  100
mat_uw <- make_pd_weight(mat_input, pn_eps = 1e-4)

Rounding the result to one decimal digit leads to

round(mat_uw, digits = 1)
#>       [,1]  [,2]  [,3]  [,4]  [,5]
#> [1,] 103.2  90.8  79.5  44.5  37.1
#> [2,]  90.8 106.5  94.2  74.1  44.5
#> [3,]  79.5  94.2 102.3  94.2  79.5
#> [4,]  44.5  74.1  94.2 106.5  90.8
#> [5,]  37.1  44.5  79.5  90.8 103.2

In case there are information about the estimation error of the different variance components, these can be specified by a weight matrix which can be specified using

mat_weight <- matrix(data = c(1000,500,20,50,200,
                               500, 1000,500,5,50,
                               20, 500, 1000,20,20,
                               50, 5, 20, 1000,200,
                               200, 50, 20, 200,1000), ncol = 5)
mat_w <- make_pd_weight(mat_input, pn_eps = 1e-4, pmat_weight = mat_weight)
round(mat_w, digits = 1)
#>       [,1]  [,2]  [,3]  [,4]  [,5]
#> [1,] 100.2  94.5  82.9  43.6  39.2
#> [2,]  94.5 100.6  94.0  60.0  45.9
#> [3,]  82.9  94.0 100.7  84.9  73.1
#> [4,]  43.6  60.0  84.9 100.3  94.2
#> [5,]  39.2  45.9  73.1  94.2 100.2

Alternatives

The Matrix package contains the function nearPD which according to its helpfile computes the nearest positive definite matrix. This function allows the specification of many more parameters. But from our experiences this can lead to the large change in a single matrix element. Furthermore, it is unclear whether it is possible to change more than one eigenvalue.

Matrix::nearPD(mat_npd)
#> $mat
#> 4 x 4 Matrix of class "dpoMatrix"
#>            [,1]      [,2]      [,3]     [,4]
#> [1,] 103.614269 75.500255 18.377845 5.013362
#> [2,]  75.500255 55.602154 12.019573 3.228358
#> [3,]  18.377845 12.019573  6.728055 1.442822
#> [4,]   5.013362  3.228358  1.442822 1.269336
#> 
#> $eigenvalues
#> [1] 1.621627e+02 4.133902e+00 9.171925e-01 1.621627e-06
#> 
#> $corr
#> [1] FALSE
#> 
#> $normF
#> [1] 10.21382
#> 
#> $iterations
#> [1] 2
#> 
#> $rel.tol
#> [1] 0
#> 
#> $converged
#> [1] TRUE
#> 
#> attr(,"class")
#> [1] "nearPD"

Session Info

sessionInfo()
#> R version 4.0.4 (2021-02-15)
#> Platform: x86_64-apple-darwin17.0 (64-bit)
#> Running under: macOS Catalina 10.15.7
#> 
#> Matrix products: default
#> BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
#> 
#> locale:
#> [1] de_CH.UTF-8/de_CH.UTF-8/de_CH.UTF-8/C/de_CH.UTF-8/de_CH.UTF-8
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] rvcetools_0.0.9
#> 
#> loaded via a namespace (and not attached):
#>  [1] bslib_0.3.0       compiler_4.0.4    pillar_1.6.3      jquerylib_0.1.4  
#>  [5] tools_4.0.4       digest_0.6.28     lattice_0.20-45   jsonlite_1.7.2   
#>  [9] evaluate_0.14     memoise_2.0.0     lifecycle_1.0.1   tibble_3.1.4     
#> [13] gtable_0.3.0      pkgconfig_2.0.3   rlang_0.4.11      Matrix_1.3-4     
#> [17] DBI_1.1.1         yaml_2.2.1        pkgdown_1.6.1     xfun_0.26        
#> [21] fastmap_1.1.0     dplyr_1.0.7       stringr_1.4.0     knitr_1.36       
#> [25] generics_0.1.0    desc_1.4.0        fs_1.5.0          sass_0.4.0       
#> [29] vctrs_0.3.8       systemfonts_1.0.2 tidyselect_1.1.1  rprojroot_2.0.2  
#> [33] grid_4.0.4        glue_1.4.2        R6_2.5.1          textshaping_0.3.5
#> [37] fansi_0.5.0       rmarkdown_2.11    tidyr_1.1.4       purrr_0.3.4      
#> [41] ggplot2_3.3.5     magrittr_2.0.1    scales_1.1.1      htmltools_0.5.2  
#> [45] ellipsis_0.3.2    assertthat_0.2.1  colorspace_2.0-2  ragg_1.1.3       
#> [49] utf8_1.2.2        stringi_1.7.4     munsell_0.5.0     cachem_1.0.6     
#> [53] crayon_1.4.1

Latest Update

2021-12-21 17:25:12 (skn)