A function to compute the biweight mean vector and covariance matrix

Compute a multivariate location and scale estimate based on Tukey's biweight weight function.

Usage

biwt.est(x, r = 0.2, med.init = robustbase::covMcd(x))

Arguments

x

a 2 x n matrix or data frame (n is the number of measurements)

r

breakdown (k/n where k is the largest number of measurements that can be replaced with arbitrarily large values while keeping the estimates bounded). Default is r = 0.2.

med.init

a (robust) initial estimate of the center and shape of the data. The format is a list with components center and cov (as in the output of covMcd() from the rrcov package). Default is the minimum covariance determinant (MCD) on the data.

Value

A list with components

biwt.mu

the final estimate of center

biwt.sig

the final estimate of shape

References

Hardin, J., Mitani, A., Hicks, L., VanKoten, B.; A Robust Measure of Correlation Between Two Genes on a Microarray, BMC Bioinformatics, 8:220; 2007.

Examples

# note that biwt.est() takes data that is 2xn where the
# goal is to find the correlation between the 2 items

samp.data <- t(MASS::mvrnorm(30,mu=c(0,0),
               Sigma=matrix(c(1,.75,.75,1),ncol=2)))

samp.bw <- biwt.est(samp.data)
samp.bw
#> $biwt.mu
#> [1] 0.2923344 0.2346656
#> 
#> $biwt.sig
#>           [,1]      [,2]
#> [1,] 0.6248218 0.5618640
#> [2,] 0.5618640 0.8092572
#> 

samp.bw.var1 <- samp.bw$biwt.sig[1,1]
samp.bw.var2 <- samp.bw$biwt.sig[2,2]
samp.bw.cov <- samp.bw$biwt.sig[1,2]

samp.bw.cor <- samp.bw.cov / sqrt(samp.bw.var1 * samp.bw.var2)
samp.bw.cor
#> [1] 0.7901505

# or:

samp.bw.cor <- samp.bw$biwt.sig[1,2] / sqrt(samp.bw$biwt.sig[1,1]*samp.bw$biwt.sig[2,2])
samp.bw.cor
#> [1] 0.7901505

##############
# to speed up the calculations, use the median/mad for the initialization:
##############

samp.init <- list()
samp.init$cov <- diag(apply(samp.data, 1, stats::mad, na.rm=TRUE))
samp.init$center <- apply(samp.data, 1, median, na.rm=TRUE)
samp.init
#> $cov
#>          [,1]     [,2]
#> [1,] 1.033435 0.000000
#> [2,] 0.000000 1.395451
#> 
#> $center
#> [1] 0.38901627 0.08568794
#> 

samp.bw <- biwt.est(samp.data, med.init = samp.init)
samp.bw.cor <- samp.bw$biwt.sig[1,2] / sqrt(samp.bw$biwt.sig[1,1]*samp.bw$biwt.sig[2,2])
samp.bw.cor
#> [1] 0.7901505