[R] Hotelling Test

[Dimitris Rizopoulos]

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
>