some time ago I've written a function for the Hotelling test, maybe
you could give it a try:
hotel.test <- function(x, y = NULL, mu = 0) {
if(!is.numeric(x) || !is.matrix(x))
stop("'x' must be a numeric matrix")
n <- nrow(x)
p <- ncol(x)
xbar <- colMeans(x, na.rm = TRUE)
if(!is.numeric(mu) || ((lmu <- length(mu)) > 1 & lmu != p))
stop("'mu' must be a numeric vector of length ", p)
if(lmu == 1)
mu <- rep(mu, p)
xbar.mu <- xbar - mu
V <- var(x, na.rm = TRUE)
out <- if(is.null(y)){
k <- n / (n - 1) * (n - p) / p
stat <- k * crossprod(xbar.mu, solve(V, xbar.mu))[1, ]
pvalue <- 1 - pf(stat, p, n - p)
list(statistic = stat, parameter = c(p, n - p), p.value =
pvalue, estimate = xbar,
null.value = mu, alternative = "two-sided", method =
"Hotelling one sample test",
data.name = deparse(substitute(x)))
} else{
if(!is.numeric(y) || !is.matrix(y))
stop("'y' must be a numeric matrix")
if(ncol(y) != p)
stop("incompatible arguments")
ny <- nrow(y)
k <- n + ny - p - 1
k1 <- (n * ny) / (n + ny)
k2 <- (n + ny - 2) * p
dif.means <- xbar - colMeans(y, na.rm = TRUE)
Vy <- var(y, na.rm = TRUE)
V <- ((n - 1) * V + (ny - 1) * Vy) / (n + ny - 2)
stat <- k * k1 * crossprod(dif.means, solve(V,
dif.means))[1, ] / k2
pvalue <- 1 - pf(stat, p, k)
list(statistic = stat, parameter = c(p, k), p.value = pvalue,
estimate = rbind(xbar, ybar = xbar - dif.means),
null.value = NULL, alternative = "two-sided", method =
"Hotelling two sample test",
data.name = c(deparse(substitute(x)),
deparse(substitute(y))))
}
class(out) <- "hotel"
out
}
print.hotel <- function(x, digits = 3, ...){
cat("\n\t", x$method, "\n\n")
if(length(dnam <- x$data.name) == 1)
cat("data:", dnam, "\n")
else
cat("data:", dnam[1], "and", dnam[2],
"\n")
pval <- if(x$p.value < 1e-04) "< 1e-04" else
paste("= ",
round(x$p.value, 4), sep = "")
cat("t = ", x$statistic, ", df1 = ", x$parameter[1],
", df2 = ",
x$parameter[2],
", p-value ", pval, "\n", sep = "")
if(!is.null(null <- x$null.value))
cat("alternative hypothesis: true mean vector is not equal
to",
paste("(", paste(round(null, digits), collapse =
",
"), ")'", sep = ""), "\n")
else
cat("alternative hypothesis: true difference in mean vectors
is not equal to 0\n")
cat("sample estimates:")
if(!is.matrix(est <- x$estimate))
cat("\nmean x-vector", paste("(", paste(round(est,
digits),
collapse = ", "), ")'", sep = ""),
"\n")
else{
rownames(est) <- c("mean x-vector", "mean
y-vector")
if(is.null(colnames(est)))
colnames(est) <- rep("", ncol(est))
print(round(est, digits))
}
cat("\n")
invisible(x)
}
# Examples
mat <- matrix(rnorm(100 * 3), 100, 3)
mat2 <- matrix(rnorm(100 * 3), 100, 3)
hotel.test(mat)
hotel.test(mat, mu = -1:1)
hotel.test(mat, y = mat2)
I hope it helps.
Best,
Dimitris
----
Dimitris Rizopoulos
Ph.D. Student
Biostatistical Centre
School of Public Health
Catholic University of Leuven
Address: Kapucijnenvoer 35, Leuven, Belgium
Tel: +32/16/336899
Fax: +32/16/337015
Web: http://www.med.kuleuven.be/biostat/
http://www.student.kuleuven.be/~m0390867/dimitris.htm
----- Original Message -----
From: "Bill Donner" <bdonner2 at yahoo.com>
To: <R-help at stat.math.ethz.ch>
Sent: Wednesday, September 07, 2005 3:48 PM
Subject: [R] Hotelling Test
> Hello R-users,
>
> I've been looking for a function performing one and two sample
> Hotelling
> test for testing equality of mean vectors. Has anyone implemented
> such a
> function in R?
>
>
> thanks a lot,
>
> Bill
>
> =============> Bill Donner
> Statistician
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
> http://www.R-project.org/posting-guide.html
>