- Se cargan los datos
D=read.table("Dados.txt",head=TRUE)
D$Trat=as.factor(D$Trat)
D$Bloco=as.factor(D$Bloco)
print(D)
## Trat Bloco X1 X2 X3 X4
## 1 T1 1 50.2 20.5 3.9 104.9
## 2 T2 1 41.8 19.5 3.7 88.0
## 3 T3 1 39.2 19.0 4.5 80.0
## 4 T4 1 33.8 20.0 4.3 80.8
## 5 T5 1 35.6 20.0 4.1 60.0
## 6 T6 1 53.4 19.2 4.2 96.4
## 7 T7 1 43.8 19.5 4.3 91.5
## 8 T8 1 50.6 19.7 4.2 91.8
## 9 T1 2 41.4 20.6 4.0 84.3
## 10 T2 2 47.2 20.1 3.6 106.5
## 11 T3 2 37.6 18.5 4.6 71.3
## 12 T4 2 49.6 20.3 4.4 106.5
## 13 T5 2 31.4 20.8 4.0 52.5
## 14 T6 2 50.2 19.5 4.5 98.8
## 15 T7 2 46.8 20.4 4.3 99.7
## 16 T8 2 57.8 19.8 4.0 84.8
## 17 T1 3 36.2 20.5 3.8 77.0
## 18 T2 3 39.6 19.3 3.6 89.8
## 19 T3 3 38.8 18.1 4.6 77.5
## 20 T4 3 35.4 20.6 4.2 83.3
## 21 T5 3 33.2 20.3 4.2 53.0
## 22 T6 3 49.6 20.3 4.3 99.1
## 23 T7 3 41.4 20.7 4.2 83.3
## 24 T8 3 41.8 20.1 4.3 70.6
## 25 T1 4 39.8 19.6 3.9 76.5
## 26 T2 4 46.6 20.1 3.7 108.7
## 27 T3 4 33.6 19.3 4.7 69.5
## 28 T4 4 41.8 20.3 4.3 95.9
## 29 T5 4 29.8 19.9 4.1 51.0
## 30 T6 4 57.8 19.9 4.5 107.2
## 31 T7 4 43.6 20.3 4.3 89.5
## 32 T8 4 46.8 20.5 4.1 81.5
- Se corre el MANOVA
m=aov(cbind(X1,X2,X3,X4)~Trat+Bloco,data=D)
anova(m)
## Analysis of Variance Table
##
## Df Pillai approx F num Df den Df Pr(>F)
## (Intercept) 1 0.99988 38097 4 18 <2e-16 ***
## Trat 7 3.12170 11 28 84 <2e-16 ***
## Bloco 3 0.54976 1 12 60 0.3602
## Residuals 21
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
- Obtención de la matriz de covarianza
SQP=crossprod(m$residuals)
CRE=SQP/m$df.residual
CRE
## X1 X2 X3 X4
## X1 18.74505952 0.016428571 -0.014702381 30.68357143
## X2 0.01642857 0.187976190 -0.003690476 0.85071429
## X3 -0.01470238 -0.003690476 0.009583333 0.09083333
## X4 30.68357143 0.850714286 0.090833333 71.79547619
- Obtención de matriz de medias
Med=aggregate(D[,3:6],by=list(D$Trat),mean)[,-1]
Med2=as.matrix(Med)
Med2
## X1 X2 X3 X4
## [1,] 41.90 20.300 3.900 85.675
## [2,] 43.80 19.750 3.650 98.250
## [3,] 37.30 18.725 4.600 74.575
## [4,] 40.15 20.300 4.300 91.625
## [5,] 32.50 20.250 4.100 54.125
## [6,] 52.75 19.725 4.375 100.375
## [7,] 43.90 20.225 4.275 91.000
## [8,] 49.25 20.025 4.150 82.175
- Obtención de las distancias entre tratamientos.
V=Med2[1,]-Med2[2,]
t(V)%*%solve(CRE)%*%V
## [,1]
## [1,] 23.68515
- Existe una libreria para realizar todo de forma más rápida
library(biotools)
## Loading required package: MASS
## ---
## biotools version 4.2
library(MASS)
Dist=D2.dist(Med2,CRE)
Dist
## 1 2 3 4 5 6 7
## 2 23.685155
## 3 63.031223 125.993410
## 4 17.571079 57.047785 26.304830
## 5 30.073588 101.542843 47.385974 37.424790
## 6 34.597979 91.840067 21.936933 17.456090 40.065822
## 7 15.533134 64.168870 23.678284 3.141215 26.009345 7.645087
## 8 29.399815 97.579518 43.372241 34.617702 15.680076 13.462747 17.252613
- A partir de la información se pueden construir dendogramas.
dendo=hclust(Dist)
plot(dendo)
# Datos de ejemplo en forma de data frame
data <- data.frame(
x1 = c(2, 3, 5, 7, 1),
x2 = c(9, 4, 6, 8, 10),
x3 = c(3, 6, 8, 1, 5)
)
# Punto de referencia (puede ser cualquier punto en el espacio)
reference_point <- c(3, 6, 8)
# Calcular la matriz de covarianza
cov_matrix <- cov(data)
# Calcular la matriz de covarianza inversa
inv_cov_matrix <- solve(cov_matrix)
# Calcular la diferencia entre cada punto y el punto de referencia
diff <- t(t(data) - reference_point)
# Calcular la distancia de Mahalanobis para cada punto
mahalanobis_distances <- sqrt(rowSums((diff %*% inv_cov_matrix) * diff))
# Mostrar las distancias de Mahalanobis
print(mahalanobis_distances)
## [1] 2.078113 1.937022 1.016288 2.784251 1.783430