Nothing
#[export]
cor_test <- function(y, x, type = "pearson", rho = 0, a = 0.05) {
n <- length(y)
if ( type == "pearson" ) {
r <- as.vector( cor(y, x) )
zh0 <- 0.5 * log( (1 + rho)/(1 - rho) )
zh1 <- 0.5 * log( (1 + r)/(1 - r) )
se <- 1/sqrt(n - 3)
}
else if ( type == "spearman" ) {
r <- as.vector( cor(Rfast::Rank(y), Rfast::Rank(x)) )
zh0 <- 0.5 * log( (1 + rho)/(1 - rho) )
zh1 <- 0.5 * log( (1 + r)/(1 - r) )
se <- 1.029563/sqrt(n - 3)
}
test <- (zh1 - zh0)/se
pvalue <- 2 * pt(abs(test), n - 3, lower.tail = FALSE)
b1 <- zh1 - qt(1 - a/2, n - 3) * se
b2 <- zh1 + qt(1 - a/2, n - 3) * se
ca <- cbind(b1, b2)
ela <- exp(2 * ca)
ci <- (ela - 1)/(ela + 1)
res <- as.vector( c(r, ci, test, pvalue) )
names(res) <- c("correlation", paste(c(a/2 * 100,
(1 - a/2) * 100), "%", sep = ""), "t-stat", "p-value")
res
}
#[export]
covar <- function(y, x) {
denom <- dim(x)[1] - 1
( Rfast::eachcol.apply(x, y) - Rfast::colmeans(x) * sum(y) ) / denom
}
#[export]
pooled.colVars <- function (x, ina, std = FALSE) {
m <- rowsum(x, ina)
m2 <- rowsum(x^2, ina)
ni <- tabulate(ina)
ni <- ni[ni > 0]
s <- (m2 - m^2/ni)
s <- Rfast::colsums(s) / (sum(ni) - length(ni) )
if (std) s <- sqrt(s)
s
}
#[export]
covlikel <- function(x, ina, a = 0.05) {
## x is the data set
## ina is a numeric vector indicating the groups of the data set
## a is the level of significance, set to 0.05 by default
ina <- as.numeric(ina)
p <- dim(x)[2] ## dimension of the data set
n <- dim(x)[1] ## total sample size
k <- max(ina) ## number of groups
nu <- tabulate(ina) ## the sample size of each group
t1 <- rep( (nu - 1)/nu, each = p^2 )
t2 <- rep(nu - 1, each = p^2 )
s <- array( dim = c(p, p, k) )
## the next 3 lines create the pooled covariance matrix
## and calculate the covariance matrix of each group
for (i in 1:k) s[, , i] <- cov( x[ina == i, ] )
mat <- t1 * s
mat1 <- t2 * s
Sp <- colSums( aperm(mat1) ) / n
deta <- apply(mat, 3, det)
pame <- det(Sp) / deta
test <- sum(nu * log(pame)) ## test statistic
dof <- 0.5 * p * (p + 1) * (k - 1) ## degrees of freedom of the asymptotic chi-square
pvalue <- pchisq(test, dof, lower.tail = FALSE) ## p-value of the test statistic
crit <- qchisq(1 - a, dof) ## critical value of the chi-square distribution
res <- c(test, pvalue, dof, crit)
names(res) <- c('test', 'p-value', 'df', 'critical')
res
}
#[export]
covmtest <- function(x, ina, a = 0.05) {
## x is the data set
## ina is a numeric vector indicating the groups of the data set
## a is the level of significance, set to 0.05 by default
p <- dim(x)[2] ## dimension of the data set
n <- dim(x)[1] ## total sample size
ina <- as.numeric(ina)
k <- max(ina) ## number of groups
nu <- tabulate(ina) ## the sample size of each group
ni <- rep(nu - 1, each = p^2)
mat <- array(dim = c(p, p, k))
## next is the covariance of each group
for (i in 1:k) mat[, , i] <- cov(x[ina == i, ])
mat1 <- ni * mat
pame <- apply(mat, 3, det) ## the detemirnant of each covariance matrix
## the next 2 lines calculate the pooled covariance matrix
Sp <- colSums( aperm(mat1) ) / ( n - k )
pamela <- det(Sp) ## determinant of the pooled covariance matrix
test1 <- sum( (nu - 1) * log(pamela/pame) )
gama1 <- ( 2 * (p^2) + 3 * p - 1 ) / ( 6 * (p + 1) * (k - 1) )
gama2 <- sum( 1/(nu - 1) ) - 1/(n - k)
gama <- 1 - gama1 * gama2
test <- gama * test1 ## this is the M (test statistic)
dof <- 0.5 * p * (p + 1) * (k - 1) ## degrees of freedom of
## the chi-square distribution
pvalue <- pchisq(test, dof, lower.tail = FALSE) ## p-value of the test statistic
crit <- qchisq(1 - a, dof) ## critical value of the chi-square distribution
result <- c(test, pvalue, dof, crit)
names(result) <- c('M.test', 'p-value', 'df', 'critical')
result
}
#[export]
covequal <- function(x, sigma, a = 0.05) {
## x is the data set
## Sigma is the assumed covariance matrix
## a is the level of significance set by default to 0.05
p <- dim(x)[2] ## dimensionality of the data
n <- dim(x)[1] ## sample size
s <- Rfast::cova(x) ## sample covariance matrix
mesa <- solve(sigma, s)
test <- n * sum( diag(mesa) ) - n * log( det(mesa) ) - n * p ## test statistic
dof <- 0.5 * p * (p + 1) ## the degrees of freedom of the chi-square distribution
pvalue <- pchisq(test, dof, lower.tail = FALSE) ## p-value of the test statistic
crit <- qchisq(1 - a, dof) ## critical value of the chi-square distribution
res <- c(test, pvalue, dof, crit)
names(res) <- c("test", "p-value", "df", "critical")
res
}
#[export]
covdist <- function(s1, s2) {
## A and B are two covariance matrices
## the order is irrelevant, s1, s2 or s2, s1 is the same
S <- solve(s1, s2)
sqrt( sum( log( eigen(S)$values )^2 ) )
}
#[export]
el.cor.test <- function(y, x, rho, tol = 1e-07) {
funa <- function(zrho, n) {
lam1 <- 0
f <- mean(zrho) - rho
der <- sum( zrho^2 / (1 + lam1 * zrho)^2 )
lam2 <- lam1 + f/der
i <- 2
while (abs(lam1 - lam2) > tol) {
i <- i + 1
lam1 <- lam2
frac <- zrho / (1 + lam1 * zrho)
f <- sum( frac ) - rho
der <- sum( frac^2 )
lam2 <- lam1 + f/der
}
list(iters = i, lam2 = lam2, p = 1 / ( 1 + lam1 * zrho ) / n )
}
x <- ( x - mean(x) ) / Rfast::Var(x, std = TRUE)
y <- ( y - mean(y) ) / Rfast::Var(y, std = TRUE)
z <- x * y
n <- length(z)
zrho <- z - rho
res <- try( funa(zrho, n), silent = TRUE )
if ( identical(class(res), "try-error") ) {
p <- iters <- NULL
info <- c(0, 10^5, 0)
} else {
p <- res$p
stat <- -2 * sum( log(n * p) )
pvalue <- pchisq(stat, 1, lower.tail = FALSE)
info <- c(res$lam2, stat, pvalue)
iters <- res$iters
}
names(info) <- c("lambda", "statistic", "p-value")
list(iters = iters, info = info, p = p)
}
#[export]
eel.cor.test <- function(y, x, rho, tol = 1e-07) {
funa <- function(zrho, zhro2, n) {
lam1 <- 0
f <- mean(zrho) - rho
der <- ( n * sum(zrho2) - sum(zrho)^2 ) / n^2
lam2 <- lam1 - f / der
i <- 2
while (abs(lam1 - lam2) > tol) {
i <- i + 1
lam1 <- lam2
com <- exp( lam1 * zrho )
up <- sum( zrho * com )
down <- sum( com )
f <- up / down - rho
der <- ( sum(zrho2 * com) * sum(com) - up^2 ) / down^2
lam2 <- lam1 - f / der
}
list(iters = i, lam2 = lam2, p = com / down )
}
x <- ( x - mean(x) ) / Rfast::Var(x, std = TRUE)
y <- ( y - mean(y) ) / Rfast::Var(y, std = TRUE)
z <- x * y
n <- length(z)
zrho <- z - rho ; zrho2 <- zrho^2
res <- try( funa(zrho, zrho2, n), silent = TRUE )
if ( identical(class(res), "try-error") ) {
p <- iters <- NULL
info <- c(0, 10^5, 0)
} else {
p <- res$p
stat <- -2 * sum( log(n * p) )
pvalue <- pchisq(stat, 1, lower.tail = FALSE)
info <- c(res$lam2, stat, pvalue)
iters <- res$iters
}
names(info) <- c("lambda", "statistic", "p-value")
list(iters = iters, info = info, p = p)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.