Astat <- function (x, tst="Wilks", univariate=TRUE, ...) 
{
    #test <- x$test
    result <- list()
    if(sum(class(x) == "anova") > 0 ) {
       result$unv <- x
       
       fss <- x[1:(nrow(x)-1),1]; ess <- x[nrow(x),1]
        peta2s <- fss/(fss+ess); f2s <- (fss/ess)
        efs <- data.frame(cbind(peta2s, f2s))
        rownames(efs) <- rownames(x)[1:(nrow(x)-1)]; colnames(efs) <- c("p.eta2", "f2s")
       result$efs <- efs
       
       return(result)
       }
    test <- tst
    repeated <- x$repeated
    ntests <- length(x$terms)
    tests <- matrix(NA, ntests, 4)
    if (!repeated) 
        SSPE.qr <- qr(x$SSPE)
    for (term in 1:ntests) {
        eigs <- Re(eigen(qr.coef(if (repeated) qr(x$SSPE[[term]]) else SSPE.qr, 
            x$SSP[[term]]), symmetric = FALSE)$values)
        tests[term, 1:4] <- switch(test, Pillai = stats:::Pillai(eigs, 
            x$df[term], x$error.df), Wilks = stats:::Wilks(eigs, 
            x$df[term], x$error.df), `Hotelling-Lawley` = stats:::HL(eigs, 
            x$df[term], x$error.df), Roy = stats:::Roy(eigs, 
            x$df[term], x$error.df))
    }
    ok <- tests[, 2] >= 0 & tests[, 3] > 0 & tests[, 4] > 0
    ok <- !is.na(ok) & ok
    tests <- cbind(x$df, tests, pf(tests[ok, 2], tests[ok, 3], 
        tests[ok, 4], lower.tail = FALSE))
    rownames(tests) <- x$terms
    colnames(tests) <- c("Df", "test stat", "approx F", "num Df", 
        "den Df", "Pr(>F)")
    tests <- structure(as.data.frame(tests), heading = paste("\nType ", 
        x$type, if (repeated) 
            " Repeated Measures", " MANOVA Tests: ", test, " test statistic\n", 
        sep = ""), class = c("anova", "data.frame"))
    result$mlt <- tests
        GG <- function(SSPE, P) {
        p <- nrow(SSPE)
        if (p < 2) 
            return(NA)
        lambda <- eigen(SSPE %*% solve(t(P) %*% P))$values
        lambda <- lambda[lambda > 0]
        ((sum(lambda)/p)^2)/(sum(lambda^2)/p)
    }
    HF <- function(gg, error.df, p) {
        ((error.df + 1) * p * gg - 2)/(p * (error.df - p * gg))
    }
    mauchly <- function(SSD, P, df) {
        if (nrow(SSD) < 2) 
            return(c(NA, NA))
        Tr <- function(X) sum(diag(X))
        p <- nrow(P)
        I <- diag(p)
        Psi <- t(P) %*% I %*% P
        B <- SSD
        pp <- nrow(SSD)
        U <- solve(Psi, B)
        n <- df
        logW <- log(det(U)) - pp * log(Tr(U/pp))
        rho <- 1 - (2 * pp^2 + pp + 2)/(6 * pp * n)
        w2 <- (pp + 2) * (pp - 1) * (pp - 2) * (2 * pp^3 + 6 * 
            pp^2 + 3 * p + 2)/(288 * (n * pp * rho)^2)
        z <- -n * rho * logW
        f <- pp * (pp + 1)/2 - 1
        Pr1 <- pchisq(z, f, lower.tail = FALSE)
        Pr2 <- pchisq(z, f + 4, lower.tail = FALSE)
        pval <- Pr1 + w2 * (Pr2 - Pr1)
        c(statistic = c(W = exp(logW)), p.value = pval)
    }
    #if (missing(test.statistic)) 
        test.statistic <- c("Pillai", "Wilks", "Hotelling-Lawley", 
            "Roy")
    test.statistic <- match.arg(test.statistic, c("Pillai", "Wilks", 
        "Hotelling-Lawley", "Roy"), several.ok = TRUE)
        
        nterms <- length(x$terms)
    if (x$repeated && univariate) {
        error.df <- x$error.df
        table <- matrix(0, nterms, 6)
        table2 <- matrix(0, nterms, 4)
        table3 <- matrix(0, nterms, 2)
        rownames(table3) <- rownames(table2) <- rownames(table) <- x$terms
        colnames(table) <- c("SS", "num Df", "Error SS", "den Df", 
            "F", "Pr(>F)")
        colnames(table2) <- c("GG eps", "Pr(>F[GG])", "HF eps", 
            "Pr(>F[HF])")
        colnames(table3) <- c("Test statistic", "p-value")
        for (term in 1:nterms) {
            SSP <- x$SSP[[term]]
            SSPE <- x$SSPE[[term]]
            P <- x$P[[term]]
            p <- ncol(P)
            PtPinv <- solve(t(P) %*% P)
            gg <- GG(SSPE, P)
            table[term, "SS"] <- sum(diag(SSP %*% PtPinv))
            table[term, "Error SS"] <- sum(diag(SSPE %*% PtPinv))
            table[term, "num Df"] <- x$df[term] * p
            table[term, "den Df"] <- error.df * p
            table[term, "F"] <- (table[term, "SS"]/table[term, 
                "num Df"])/(table[term, "Error SS"]/table[term, 
                "den Df"])
            table[term, "Pr(>F)"] <- pf(table[term, "F"], table[term, 
                "num Df"], table[term, "den Df"], lower.tail = FALSE)
            table2[term, "GG eps"] <- gg
            table2[term, "HF eps"] <- HF(gg, error.df, p)
            table3[term, ] <- mauchly(SSPE, P, x$error.df)
        }
        attr(table, "heading") <- paste("\nUnivariate Type", x$type, "Repeated-Measures ANOVA Assuming Sphericity\n")
        class(table)<- c("anova", "array")
        
        fss <- table[,1]; ess <- table[,3]
        peta2s <- fss/(fss+ess); f2s <- (fss/ess)
        efs <- data.frame(cbind(peta2s, f2s))
        rownames(efs) <- rownames(table); colnames(efs) <- c("p.eta2", "f2s")
        
        result$unv <- table        
        result$efs <- efs
        
        table3 <- na.omit(table3)
        if (nrow(table3) > 0) {
            attr(table3, "heading") <- paste("\nMauchly Tests for Sphericity\n")
            class(table3)<- c("anova", "array")
            result$mct <- table3
            #print.anova(table3)
            #cat("\n\nGreenhouse-Geisser and Huynh-Feldt Corrections\n", 
            #    "for Departure from Sphericity\n\n")
            table2[, "Pr(>F[GG])"] <- pf(table[, "F"], table2[, 
                "GG eps"] * table[, "num Df"], table2[, "GG eps"] * 
                table[, "den Df"], lower.tail = FALSE)
            table2[, "Pr(>F[HF])"] <- pf(table[, "F"], pmin(1, 
                table2[, "HF eps"]) * table[, "num Df"], pmin(1, 
                table2[, "HF eps"]) * table[, "den Df"], lower.tail = FALSE)
            table2 <- na.omit(table2)
            table2.12 <- table2[, 1:2, drop = FALSE]
            class(table2.12)<- c("anova", "array")
            attr(table2.12, "heading") <- "\nGreenhouse-Geisser corrections"
            result$gg <- table2.12
            #print.anova(table2[, 1:2, drop = FALSE])
            cat("\n")
            table2.34 <- table2[, 3:4, drop = FALSE]
            class(table2.34)<- c("anova", "array")
            attr(table2.34, "heading") <- "\nHuynh-Feldt corrections"
            result$hf <- table2.34
            #print.anova(table2[, 3:4, drop = FALSE])
            if (any(table2[, "HF eps"] > 1)) 
                warning("HF eps > 1 treated as 1")
        }
        }
    return(result)
    #print(tests)
    #invisible(x)
}