Métodos multivariados
Antonio Ortiz y Dylan Guevara
26/6/2020
#a
library("plot3D")
xt=c(5, 1, 3)
yt=c(-1, 3, 1)
m <- cbind(xt, yt)
arrows3D(x0 = c(0,0), y0 = c(0,0), z0 = c(0,0), x1 =m[1,], y1 = m[2,], z1 = m[3,], col = c("green", "blue"), lwd = 2, d = 3, main = "Vectors 3D", bty ="g", ticktype = "detailed")
text3D(m[1,], m[2,], m[3,], c("xt", "yt"), col = c("black", "black"), add=TRUE, colkey = FALSE)
#b
x <- t(xt)
y <- t(yt)
l <- sqrt(sum(x*x)); norm(x, type = "2")## [1] 5.91608l <- sqrt(sum(y*y)); norm(x, type = "2")## [1] 5.91608angulo = function(u, v){
prod.punto=t(u)%*%v; prod.punto
norma.u <- norm(u, type = "2"); norma.u
norma.v <- norm(v, type = "2"); norma.v
phi <- acos(prod.punto / (norma.u * norma.v)); phi
}
angulo(xt,yt)## [,1]
## [1,] 1.51981#c
x1 <- x-3
y1 <- y-1
m1 <- rbind(x1, y1)
arrows3D(x0 = c(0,0), y0 = c(0,0), z0 = c(0,0), x1 =m1[,1], y1 = m1[,2], z1 = m1[,3], col = c("green", "blue"), lwd = 2, d = 3, main = "Vectors 3D", bty ="g", ticktype = "detailed")
text3D(m[1,], m[2,], m[3,], c("x1", "y1"), col = c("black", "black"), add=TRUE, colkey = FALSE)
Ejercicio 1: Un productor decide probar el funcionamiento de su máquina y para ello, luego de cosechar una parcela, cuenta en 10 unidades de 1 m2 la cantidad de semillas que quedan en el suelo. Las normas técnicas indican que la media del número de semillas caídas por m2 no debería ser superior a 80.
- Construir un intervalo de confianza para u con una confianza del 90%.
- Concluir sobre el funcionamiento de la máquina
\[Ho: \bar{x} = 80 sem.m^{-2}\\Ha: \bar{x} \ne 80 sem.m^{-2}\]
datos <- c(77,73,82,82,79,81,78,76,76,75)
summary(datos)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 73.0 76.0 77.5 77.9 80.5 82.0test <- t.test(x = datos, mu = 80, alternative = 'l',conf.level = 0.90)
test$p.value## [1] 0.02943033ifelse(test$p.value<0.1,
'Rechazo Ho: menos de 80 sem por m^2',
'No rechazo Ho: más de 80 sem por m^2')## [1] "Rechazo Ho: menos de 80 sem por m^2"Conclusión:El funcionamiento de la máquina es bueno. La cantidad de semillas por metro cuadrado es menor a 80.
Ejercicio 2. El espárrago es una planta perenne cuyo cultivo comercial puede tener una duración de 15 años y su implantación es costosa. Dada la extensión del sistema radicular, la profundidad del suelo es fundamental, considerándose indispensable contar con un promedio mínimo de 80 centímetros de sustrato permeable. Se realizan 14 determinaciones de la profundidad del sustrato permeable (en cm) en puntos tomados al azar en dos campos (A y B).
- A partir de los intervalos de confianza al 95% determinar si estos campos son aptos para el cultivo.
\[Ho: \bar{x} = 80 cm\\Ha: \bar{x} > 80 cm\] b) ¿Hay diferencias en la profundidad del sustrato permeable entre ambos campos?
\[Ho: \mu_{A} = \mu_{B}\\Ha:\mu_{A}\ne \mu_{B}\]
Muestra <- c(1,2,3,4,5,6,7,8,9,10,11,12,13,14)
A <- c(72,78,86,78,90,104,76,70,83,75,90,81,85,72)
B <- c(86,90,76,76,82,89,93,81,83,97,108,98,90,83)
dfAB <- data.frame(Muestra, A, B)
summary(dfAB)## Muestra A B
## Min. : 1.00 Min. : 70.00 Min. : 76.00
## 1st Qu.: 4.25 1st Qu.: 75.25 1st Qu.: 82.25
## Median : 7.50 Median : 79.50 Median : 87.50
## Mean : 7.50 Mean : 81.43 Mean : 88.00
## 3rd Qu.:10.75 3rd Qu.: 85.75 3rd Qu.: 92.25
## Max. :14.00 Max. :104.00 Max. :108.00ptA <- t.test(x = A, mu = 80, alternative = 'g',conf.level = 0.95)
ptA$p.value## [1] 0.2851637ifelse(ptA$p.value<0.05,
'Rechazo Ho: prufundidad mayor a 80 cm',
'No rechazo Ho: profundidad menor a 80 cm')## [1] "No rechazo Ho: profundidad menor a 80 cm"ptB <- t.test(x = B, mu = 80, alternative = 'g',conf.level = 0.95)
ptB$p.value## [1] 0.002643557ifelse(ptB$p.value<0.05,
'Rechazo Ho, Profundidad mayor a 80 cm',
'No rechazo Ho, profundidad menor a 80 cm')## [1] "Rechazo Ho, Profundidad mayor a 80 cm"ptAB <- t.test(x = A, y = B, mu = 0, alternative = 't',conf.level = 0.95)
ptAB$p.value## [1] 0.06622161ifelse(ptAB$p.value<0.05,
'Rechazo Ho: Las medias no son iguales',
'No rechazo Ho: medias iguales')## [1] "No rechazo Ho: medias iguales"boxplot(A,B,names=c("Campo A","Campo B"))
medias <- c(mean(A),mean(B))
points(medias,pch=18,col="red")
- Verificando varianzas
\[H_o: \sigma^2_{A} = \sigma^2_{B} \\H_a: \sigma^2_{A} \ne \sigma^2_{B}\]
pvar = var.test(x = A, y = B, ratio = 1, alternative = 't', conf.level = 0.95)
pvar$p.value## [1] 0.9295384ifelse(ptAB$p.value<0.05,
'Rechazo Ho: varianzas desiguales',
'No rechazo Ho: varianzas iguales')## [1] "No rechazo Ho: varianzas iguales"Conclusión:
- Mediante la prueba T-Student con un intervalo de confianza al 95%, se determinó que el campo A no es apto y que el campo B es apto para el cultivo de espárrago.
- De acuedo a la prueba T-Student, no hay diferencia estadística en la profundidad del sustrato permeable entre el campo A y el campo B. Sin embargo, el p-valor es de 0,066.
Prueba t-Dos muestras dependientes
Suponga que se hacen dos cortes sucesivos (separados 30 días) en el cultivo de limonaria y se extrae su aceite esencial. Para eso se diseñó un experimento con 12 parcelas cada una con 2 subparcelas bien definidas con el objeto de hacer el corte en cada subparcela sin que una interfiera en la otra, es decir, el primer corte se hizo en la primera subparcela de la misma parcela, donde a los 30 días se hizo el segundo corte.
AE1 <- c(3.12, 3.12, 3.15, 3.18, 3.22, 3.26, 3.2, 3.18)
AE2 <- c(3.11, 3.1, 3.11, 3.14, 3.15, 3.18, 3.15, 3.18)
dfAE <- data.frame(AE1, AE2); dfAEView(dfAE)
summary(dfAE)## AE1 AE2
## Min. :3.120 Min. :3.100
## 1st Qu.:3.143 1st Qu.:3.110
## Median :3.180 Median :3.145
## Mean :3.179 Mean :3.140
## 3rd Qu.:3.205 3rd Qu.:3.158
## Max. :3.260 Max. :3.180var(dfAE)## AE1 AE2
## AE1 0.002355357 0.0012714286
## AE2 0.001271429 0.0009714286boxplot(dfAE, col="white")
medias <- c(mean(AE1),mean(AE2))
points(medias,pch=18,col="red")
cor(dfAE$AE1, dfAE$AE2, method ="pearson")## [1] 0.8405395## suponiendo que son no dependientes ##
ptAE <- t.test(dfAE$AE1, dfAE$AE2, alternative = "t", mu = 0,
var.equal = T, conf.level = 0.95); ptAE##
## Two Sample t-test
##
## data: dfAE$AE1 and dfAE$AE2
## t = 1.9002, df = 14, p-value = 0.0782
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.004987255 0.082487255
## sample estimates:
## mean of x mean of y
## 3.17875 3.14000ptAE$p.value## [1] 0.07819718ifelse(ptAE$p.value<0.05,
"Rechazo Ho: Las medias no son iguales",
'No rechazo Ho: medias iguales')## [1] "No rechazo Ho: medias iguales"## prueba de varianzas mediante prueba f
pvAE <- var.test(dfAE$AE1, dfAE$AE2, ratio = 1, alternative = "t",
conf.level = 0.95)
pvAE$p.value## [1] 0.2654353ifelse(pvAE$p.value<0.05,
"Rechazo Ho: var desiguales",
'No rechazo Ho: var iguales')## [1] "No rechazo Ho: var iguales"## pareadas ##
ptAE_p <- t.test(dfAE$AE1, dfAE$AE2, paired = T, alternative = "t", mu = 0,
var.equal = T, conf.level = 0.95)
ptAE_p$p.value## [1] 0.005789683ifelse(ptAE_p$p.value<0.05,
"Rechazo Ho: Las medias no son iguales",
'No rechazo Ho: medias iguales')## [1] "Rechazo Ho: Las medias no son iguales"verde <- c(9.33, 8.74, 9.31, 8.27, 10.22, 10.13, 10.42, 10.62, 15.25, 16.22, 17.24,
12.77, 12.07,11.03,12.48,12.12,15.38,14.21,9.69,14.35,38.71,44.74,36.67,
37.21,8.73,7.94,8.37,7.86,8.45,6.79,8.34,7.54,14.04,13.51,13.33,12.77)
rojo <- c(19.14,19.55,19.24,16.37,25,25.32,27.12,26.28,38.89,36.37,40.74,67.50,
33.03,32.37,31.31,33.33,40,40.48,37.9,40.15,77.14,78.57,71.43,45,23.27,
20.87,21.16,21.78,26.32,22.73,26.77,24.87,44.44,37.93,37.93,60.87)
sp <- gl(3,12,36, labels = c("SS", "JL", "LP"))
loc <- gl(2,6,36, labels = c("L1", "L2"))
data <- data.frame(verde, rojo, sp)
Y <- cbind(verde, rojo)
modmanova <- manova(Y~sp)
summary(modmanova, test = "Wilks")## Df Wilks approx F num Df den Df Pr(>F)
## sp 2 0.6767 3.4501 4 64 0.01291 *
## Residuals 33
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1resmanova <- (modmanova$residuals)
royston::royston.test(resmanova)## $H
## [1] 25.85221
##
## $p.value
## [1] 2.115054e-06
##
## $distribution
## [1] "Data analyzed have a non-normal distribution."library(biotools)## Loading required package: rpanel## Loading required package: tcltk## Package `rpanel', version 1.1-4: type help(rpanel) for summary information## Loading required package: tkrplot## Loading required package: MASS## Loading required package: lattice## Loading required package: SpatialEpi## Loading required package: sp## ---
## biotools version 3.1## 
boxplot(rojo~sp)
boxplot(verde~sp)
#ANOVA POR VARIABLE
anovaverde <- aov(verde~sp, data=data)
summary(anovaverde)## Df Sum Sq Mean Sq F value Pr(>F)
## sp 2 965.2 482.6 7.415 0.00219 **
## Residuals 33 2147.7 65.1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1resverde <- anovaverde$residuals
shapiro.test(resverde)##
## Shapiro-Wilk normality test
##
## data: resverde
## W = 0.88653, p-value = 0.001491bartlett.test(resverde,data$sp)##
## Bartlett test of homogeneity of variances
##
## data: resverde and data$sp
## Bartlett's K-squared = 32.714, df = 2, p-value = 7.876e-08anovarojo <- aov(rojo~sp, data=data)
summary(anovarojo)## Df Sum Sq Mean Sq F value Pr(>F)
## sp 2 2125 1062.6 4.697 0.016 *
## Residuals 33 7465 226.2
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1resrojo <- anovarojo$residuals
shapiro.test(resrojo)##
## Shapiro-Wilk normality test
##
## data: resrojo
## W = 0.81736, p-value = 3.672e-05bartlett.test(resrojo,sp)##
## Bartlett test of homogeneity of variances
##
## data: resrojo and sp
## Bartlett's K-squared = 1.5991, df = 2, p-value = 0.4495TukeyHSD(anovarojo, "sp")## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = rojo ~ sp, data = data)
##
## $sp
## diff lwr upr p adj
## JL-SS 16.5991667 1.532054 31.6662795 0.0282769
## LP-SS 0.6183333 -14.448779 15.6854461 0.9944259
## LP-JL -15.9808333 -31.047946 -0.9137205 0.0357363plot(TukeyHSD(anovarojo, "sp"))
TukeyHSD(anovaverde, "sp")## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = verde ~ sp, data = data)
##
## $sp
## diff lwr upr p adj
## JL-SS 10.01167 1.930131 18.093202 0.0124725
## LP-SS -1.73750 -9.819035 6.344035 0.8584379
## LP-JL -11.74917 -19.830702 -3.667631 0.0031513plot(TukeyHSD(anovaverde, "sp"))
datab <- data.frame(verde, rojo, sp, loc)
Yb <- cbind(verde, rojo)
manovab <- manova(Yb~sp+sp*loc, data=datab)
summary(manovab)## Df Pillai approx F num Df den Df Pr(>F)
## sp 2 0.48953 4.8613 4 60 0.0018357 **
## loc 1 0.46930 12.8225 2 29 0.0001024 ***
## sp:loc 2 0.28751 2.5184 4 60 0.0504951 .
## Residuals 30
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 15. Análisis de componentes principales (ACP)
Se tomaron 40 muestras de moniliaisis del cacao (“Moniliphthora roreri”) en cultivos de Ecuador, Colombia, Perú, Venezuela y algunos países de Centroamérica. A las muestras de monilia se les realizó un genotipado con 11 marcadores moleculares tipo microsatélite. A partir de esta genotipificación se contruto una matriz de datos binarios donde 1 y 0 representan la presencia y ausencia de un alelo, respetivamente. En total se encontraron 38 alelos en las 40 muestras de moniliaisis, cada alelos se consideró como una variable para los análisis. Los datos se tomaron del artículo “The cacao pathogen Moniliophthora roreri (Marasmiaceae) possesses biallelic A and B mating loci but reproduces clonally” de Díaz y Aime (2016).
Se realizó un análisis de clúster para obtener un dendograma, el cual fue soportado por un Bootstrap con 1000 remuestreos. Luego se utilizó un análisis de componentes principales para observar la relación entre muestras.
library(readxl)
acp <- read_excel("C:/Users/Roxsana/Desktop/Métodos multivariados/acp.xlsx")
data <- data.frame(acp)
ndata <- data.frame(data[,-1], row.names = data$id)5.1. Análisis de Cluster
nndata <- ndata[-1]
matrix <- dist(nndata, method="binary")
clust <- hclust(matrix, method="complete")
plot(clust, cex=0.5, hang=-1)
#Bootstrapping
datost<- t(nndata)
library(pvclust)
boot <- pvclust(datost, method.dist="binary", method.hclust="complete", nboot=1000)## Bootstrap (r = 0.49)... Done.
## Bootstrap (r = 0.59)... Done.
## Bootstrap (r = 0.68)... Done.
## Bootstrap (r = 0.78)... Done.
## Bootstrap (r = 0.89)... Done.
## Bootstrap (r = 1.0)... Done.
## Bootstrap (r = 1.08)... Done.
## Bootstrap (r = 1.19)... Done.
## Bootstrap (r = 1.3)... Done.
## Bootstrap (r = 1.38)... Done.plot(boot, cex=0.5)
5.2. Análisis de componentes principales
library(scatterplot3d)
colour=as.character(ndata[,1])
pca<-princomp(ndata[,-1])
scatterplot3d(pca$score[,1:3], color=c(colour),
pch=19, box=0, type="h", angle=100,
grid = TRUE, ylab = "PC 2 (19%)", xlab = "PC 1 (38.4 %)", zlab = "PC3 (6.5%)")
Conclusión:
Mediante el análisis de clúster se observa la conformación de dos grandes grupos, estos dos grandes grupos se forman por la influencia del factor sexo (datos no mostrados), las muestras de moniliasiis de cada uno de los grupos tiene diferente sexo.
5.3. Análisis de conglomerados
library(factoextra)## Loading required package: ggplot2## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBaset.seed(123)
d2f=data.frame(ndata[,-c(2:3)])
km_clusters <- kmeans(x = d2f[,-1], centers = 2, nstart = 50)
cluster_mon <- km_clusters$cluster
fviz_cluster(object = km_clusters, data = d2f[,-1], show.clust.cent = TRUE,
ellipse.type = "euclid", star.plot = TRUE, repel = TRUE,
pointsize=0.5,outlier.color="darkred", panel_rows(km_clusters)) +
labs(title = "Agrupamiento K-Means") +
theme_bw() + theme(legend.position = "right")
Conclusión:
Al igual que en el análisis de clúster se observa la conformación de dos grandes grupos determinados por el sexo de los aislamientos de monilia.
library(dplyr)##
## Attaching package: 'dplyr'## The following object is masked from 'package:MASS':
##
## select## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recodedf=data.frame(iris)
iris.2 <- iris[,-5]
species <- iris[,5]
di=data.frame(df) %>%
mutate(Species=dplyr::recode(Species,
setosa="st",
versicolor="vs",
virginica="vg"))
pairs(di[,1:4], col = species,lower.panel = NULL)
par(xpd = TRUE)
legend(x = 0.05, y = 0.4, cex = 2,
legend=as.character(levels(species)),
fill = unique(species))
par(xpd = NA)
set.seed(20)
k.means.fit <- kmeans(df[,1:4], 3)
k.means.fit$centers## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 1 5.006000 3.428000 1.462000 0.246000
## 2 5.901613 2.748387 4.393548 1.433871
## 3 6.850000 3.073684 5.742105 2.071053k.means.fit$ifault## [1] 0grupos=k.means.fit$cluster
table(di$Species,grupos) #matriz de confusión (puedo hacerlo porque conozco las variedades)## grupos
## 1 2 3
## st 50 0 0
## vs 0 48 2
## vg 0 14 36dif=data.frame(di,grupos)
dif=data.frame(dif) %>%
mutate(grupos=dplyr::recode(grupos,
"3"="st",
"2"="vs",
"1"="vg"))
table(dif$grupos,dif$Species)##
## st vs vg
## st 0 2 36
## vg 50 0 0
## vs 0 48 14library(factoextra)
d2 <- scale(di[,1:4])
rownames(d2) <- di$Species
#www:for total within sum of square
fviz_nbclust(x = d2, FUNcluster = kmeans, method = "wss", k.max = 15,
diss = get_dist(d2, method = "euclidean"), nstart = 50)
set.seed(123)
d2f=data.frame(d2)
km_clusters <- kmeans(x = d2f, centers = 3, nstart = 50)
# Las funciones del paquete factoextra emplean el nombre de las filas del
# dataframe que contiene los datos como identificador de las observaciones.
# Esto permite añadir labels a los gráficos.
fviz_cluster(object = km_clusters, data = d2f, show.clust.cent = TRUE,
ellipse.type = "euclid", star.plot = TRUE, repel = TRUE,
pointsize=0.5,outlier.color="darkred") +
labs(title = "Resultados clustering K-means") +
theme_bw() + theme(legend.position = "none")
require(cluster)## Loading required package: clusterpam.res <- pam(d2f, 3)
# Visualización
fviz_cluster(pam.res, geom = "point", ellipse.type = "norm",
show.clust.cent = TRUE,star.plot = TRUE)+
labs(title = "Resultados clustering K-means")+ theme_bw()
# Load iris dataset
data(iris)
# PCA
pca <- prcomp(iris[,-5], scale=TRUE)
df.pca <- pca$x
# Cluster over the three first PCA dimensions
kc <- kmeans(df.pca[,1:3], 4)
fviz_pca_biplot(pca, label="var", habillage=as.factor(kc$cluster)) +
labs(color=NULL) + ggtitle("") +
theme(text = element_text(size = 15),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black"),
legend.key = element_rect(fill = "white"))
library(cluster)
fviz_nbclust(x = d2f[,1:4], FUNcluster = pam, method = "wss", k.max = 15,
diss = dist(d2, method = "manhattan"))
set.seed(123)
pam_clusters <- pam(x = d2f[,1:4], k = 4, metric = "manhattan")
pam_clusters$medoids## Sepal.Length Sepal.Width Petal.Length Petal.Width
## st.7 -1.0184372 0.7861738 -1.27910398 -1.3110521
## vg.12 1.1553023 -0.1315388 0.98680210 1.1816087
## vs.28 0.1891958 -0.3609670 0.42032558 0.3944526
## vs.19 -0.2938574 -1.2786796 0.08043967 -0.1303181fviz_cluster(object = pam_clusters, data = d2f[,1:4],
ellipse.type = "t",repel = TRUE) +
theme_bw() + labs(title = "Resultados clustering PAM") +
theme(legend.position = "none")
medoids <- prcomp(d2f[,1:4])$x
medoids <- medoids[rownames(pam_clusters$medoids), c("PC1", "PC2")]
medoids <- as.data.frame(medoids)
colnames(medoids) <- c("x", "y")
# Creación del gráfico
fviz_cluster(object = pam_clusters, data = d2f[,1:4], ellipse.type = "t",
repel = TRUE) + theme_bw() +
# Se resaltan las observaciones que actúan como medoids
geom_point(data = medoids, color = "firebrick", size = 2) +
labs(title = "Resultados clustering PAM") +
theme(legend.position = "none")
### método clara
library(simstudy)
library(data.table)##
## Attaching package: 'data.table'## The following objects are masked from 'package:dplyr':
##
## between, first, lastset.seed(555)
dt1 <- genCorData(200, mu = c(37.88,16.79, 11.74,21.4,34.04),
sigma = 1, rho = 0.5, corstr = "cs" )
dt2 <- genCorData(200, mu = c(35.94,17.90, 11.81,21.77,32.86),
sigma = 1, rho = 0.7, corstr = "cs" )
dac <- rbind(dt1,dt2)
DDA=gl(2,200,400,labels=c("DDA224","DDA231"))
dfc=data.frame(DDA,dac)
View(dfc)
colnames(dfc) <- c("DDA","id","L", "a","b","C","h")
dfc=dfc[,-2]
head(dfc)dac <- rbind(dt1,dt2)
DDA=gl(2,200,400,labels=c("DDA224","DDA231"))
dfc=data.frame(DDA,dac)
colnames(dfc) <- c("DDA","id","L", "a","b","C","h")
dfc=dfc[,-2] # quito id
head(dfc)options(digits = 3)
library(tidyverse)## -- Attaching packages ------------------------- tidyverse 1.3.0 --## v tibble 3.0.1 v purrr 0.3.3
## v tidyr 1.1.0 v stringr 1.4.0
## v readr 1.3.1 v forcats 0.5.0## -- Conflicts ---------------------------- tidyverse_conflicts() --
## x data.table::between() masks dplyr::between()
## x dplyr::filter() masks stats::filter()
## x data.table::first() masks dplyr::first()
## x dplyr::lag() masks stats::lag()
## x data.table::last() masks dplyr::last()
## x car::recode() masks dplyr::recode()
## x dplyr::select() masks MASS::select()
## x purrr::some() masks car::some()
## x purrr::transpose() masks data.table::transpose()dfc %>%
split(.$DDA) %>%
map(select,-c(DDA)) %>%
map(cor) ## $DDA224
## L a b C h
## L 1.000 0.468 0.415 0.535 0.417
## a 0.468 1.000 0.463 0.505 0.454
## b 0.415 0.463 1.000 0.391 0.437
## C 0.535 0.505 0.391 1.000 0.526
## h 0.417 0.454 0.437 0.526 1.000
##
## $DDA231
## L a b C h
## L 1.000 0.664 0.681 0.699 0.680
## a 0.664 1.000 0.714 0.683 0.735
## b 0.681 0.714 1.000 0.686 0.696
## C 0.699 0.683 0.686 1.000 0.707
## h 0.680 0.735 0.696 0.707 1.000library(cluster)
library(factoextra)
clara_clusters <- clara(x = dfc, k = 2, metric = "manhattan", stand = TRUE,
samples = 60, pamLike = TRUE)
table(clara_clusters$clustering)##
## 1 2
## 197 203fviz_cluster(object = clara_clusters, ellipse.type = "t", geom = "point",
pointsize = 1.5) + theme_bw() +
labs(title = "Resultados clustering CLARA") +
theme(legend.position = "none")
medoides <- prcomp(dfc[,-1])$x
medoides <- medoides[rownames(pam_clusters$medoides), c("PC1", "PC2")]
medoides <- as.data.frame(medoides)
colnames(medoides) <- c("x", "y")
# Creación del gráfico
fviz_cluster(object = clara_clusters, data = d2f[,-1], ellipse.type = "t",
repel = TRUE,pointsize = 0.5) + theme_bw() +
geom_point(data = medoides, color = "firebrick", size = 1,) +
labs(title = "Resultados clustering PAM") +
theme(legend.position = "none")
#Hacer para el método clara el BIplot, para saber cual de las variables quedo mejor representadas.
pcacl <- prcomp(dfc[1:200,2:6], scale=TRUE)
df.pca1 <- pcacl$x
clara_clusterspca <- clara(x = df.pca1, k = 2, metric = "manhattan", stand = TRUE,
samples = 60, pamLike = TRUE)
table(clara_clusters$clustering)##
## 1 2
## 197 203fviz_pca_biplot(pcacl, label="var", habillage=as.factor(clara_clusterspca$clustering)) +
labs(color=NULL) + ggtitle("") +
theme(text = element_text(size = 15),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black"),
legend.key = element_rect(fill = "white"))
#segundo muestreo
pcacl2 <- prcomp(dfc[201:400,2:6], scale=TRUE)
df.pca2 <- pcacl$x
clara_clusterspca2 <- clara(x = df.pca2, k = 2, metric = "manhattan", stand = TRUE,
samples = 60, pamLike = TRUE)
table(clara_clusterspca2$clustering)##
## 1 2
## 94 106fviz_pca_biplot(pcacl2, label="var", habillage=as.factor(clara_clusterspca2$clustering)) +
labs(color=NULL) + ggtitle("") +
theme(text = element_text(size = 15),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black"),
legend.key = element_rect(fill = "white"))