r manova花生两种方式

Posted 2021-05-24

X_df <- read.csv("peanut.csv")
X_df$Location <- factor(X_df$Location)
X_df$Variety <- factor(X_df$Variety)
X <- as.matrix(X_df[, 3:5])
## Summary statistics.
x_bar <- colMeans(X)
x_bar_l <- by(X, X_df$Location, colMeans)
x_bar_k <- by(X, X_df$Variety, colMeans)
Loc_Var <- factor(rep(c("1_1", "2_1", "1_2", "2_2", "1_3", "2_3"), each = n))
x_bar_lk <- by(X, Loc_Var, colMeans)
S_lk <- by(X, Loc_Var, var)
## MANOVA calculations.
SSP_1 <- SSP_2 <- SSP_int <- matrix(0, nrow = p, ncol = p)
for(l in 1:g)
SSP_1 <- SSP_1 + b * n * (x_bar_l[[l]] - x_bar) %*% t(x_bar_l[[l]] - x_bar)
for(k in 1:b)
SSP_2 <- SSP_2 + g * n * (x_bar_k[[k]] - x_bar) %*% t(x_bar_k[[k]] - x_bar)
for(l in 1:g) {
for(k in 1:b) {
innards <- x_bar_lk[[k + (l - 1) * b]] - x_bar_l[[l]] - x_bar_k[[k]] + x_bar
SSP_int <- SSP_int + n * innards %*% t(innards)
}
}
SSP_res <- Reduce("+", S_lk)

####################OR
for(l in 1:g) {
for(k in 1:b) {
for(r in 1:n){
SSP_res <- SSP_res + t(as.matrix(dta3[(l - 1) * 3 * n + (k - 1) * n + r, 3:5]) - x_bar_lk[(l - 1) * 3 + k, , drop = FALSE]) %*%
(as.matrix(dta3[(l - 1) * 3 * n + (k - 1) * n + r, 3:5]) - x_bar_lk[(l - 1) * 3 + k, , drop = FALSE])
SSP_cor <- SSP_cor + t(as.matrix(dta3[(l - 1) * 3 * n + (k - 1) * n + r, 3:5]) - x_bar) %*% (as.matrix(dta3[(l - 1) * 3 * n + (k - 1) * n + r, 3:5]) - x_bar)
}
}
}
######################
## Testing interaction.
Lambda_int <- det(SSP_res) / det(SSP_int + SSP_res)
Lambda_scaled <- (((g * b * (n - 1) - p + 1) / 2) /
((abs((g - 1) * (b - 1) - p) + 1) / 2)) * (1 - Lambda_int) / Lambda_int
p_val_int <- 1 - pf(Lambda_scaled, abs((g - 1) * (b - 1) - p) + 1, g * b * (n - 1) - p + 1)

## Testing location main effect.
Lambda_1 <- det(SSP_res) / det(SSP_1 + SSP_res)
Lambda_scaled <- (((g * b * (n - 1) - p + 1) / 2) / ((abs((g - 1) - p) + 1) / 2)) *
(1 - Lambda_1) / Lambda_1
p_val_1 <- 1 - pf(Lambda_scaled, abs((g - 1) - p) + 1, g * b * (n - 1) - p + 1)

## Testing variety main effect.
Lambda_2 <- det(SSP_res) / det(SSP_2 + SSP_res)
Lambda_scaled <- (((g * b * (n - 1) - p + 1) / 2) / ((abs((b - 1) - p) + 1) / 2)) *
(1 - Lambda_2) / Lambda_2
p_val_2 <- 1 - pf(Lambda_scaled, abs((b - 1) - p) + 1, g * b * (n - 1) - p + 1)
summary(manova(cbind(X_1, X_2, X_3) ~ Location * Variety, data = X_df), test = "Wilks")