En este documento se presenta el flujo de trabajo utilizado para
procesar los espectros compartidos en las carpetas 3_T 50 y
3 SST50
Primero, creamos una lista con los nombres de los archivos contenidos
en esta carpeta, para poder importarlos a R:
nombres <- list.files(, pattern = 'CSV')
nombres## [1] "@A_50C_LF1_1.CSV" "@A_50C_LF1_2.CSV" "@A_50C_LF1_3.CSV"
## [4] "@A_50C_LF2_1.CSV" "@A_50C_LF2_2.CSV" "@A_50C_LF2_3.CSV"
## [7] "@A_50C_LF3_1.CSV" "@A_50C_LF3_2.CSV" "@A_50C_LF3_3.CSV"
## [10] "@A_50C_LF4_1.CSV" "@A_50C_LF4_2.CSV" "@A_50C_LF4_3.CSV"
## [13] "@A_50C_LF5_1.CSV" "@A_50C_LF5_2.CSV" "@A_50C_LF5_3.CSV"
## [16] "@A_50C_UF1_1.CSV" "@A_50C_UF1_2.CSV" "@A_50C_UF1_3.CSV"
## [19] "@A_50C_UF2_1.CSV" "@A_50C_UF2_2.CSV" "@A_50C_UF2_3.CSV"
## [22] "@A_50C_UF3_1.CSV" "@A_50C_UF3_2.CSV" "@A_50C_UF3_3.CSV"
## [25] "@A_50C_UF4_1.CSV" "@A_50C_UF4_2.CSV" "@A_50C_UF4_3.CSV"
## [28] "@A_50C_UF5_1.CSV" "@A_50C_UF5_2.CSV" "@A_50C_UF5_3.CSV"
## [31] "@A_SS_50C_LF1_1.CSV" "@A_SS_50C_LF1_2.CSV" "@A_SS_50C_LF1_3.CSV"
## [34] "@A_SS_50C_LF2_1.CSV" "@A_SS_50C_LF2_2.CSV" "@A_SS_50C_LF2_3.CSV"
## [37] "@A_SS_50C_LF3_1.CSV" "@A_SS_50C_LF3_2.CSV" "@A_SS_50C_LF3_3.CSV"
## [40] "@A_SS_50C_LF4_1.CSV" "@A_SS_50C_LF4_2.CSV" "@A_SS_50C_LF4_3.CSV"
## [43] "@A_SS_50C_LF5_1.CSV" "@A_SS_50C_LF5_2.CSV" "@A_SS_50C_LF5_3.CSV"
## [46] "@A_SS_50C_UF1_1(1).CSV" "@A_SS_50C_UF1_1.CSV" "@A_SS_50C_UF1_2.CSV"
## [49] "@A_SS_50C_UF1_3.CSV" "@A_SS_50C_UF2_1.CSV" "@A_SS_50C_UF2_2.CSV"
## [52] "@A_SS_50C_UF2_3.CSV" "@A_SS_50C_UF3_1.CSV" "@A_SS_50C_UF3_2.CSV"
## [55] "@A_SS_50C_UF3_3.CSV" "@A_SS_50C_UF4_1.CSV" "@A_SS_50C_UF4_2.CSV"
## [58] "@A_SS_50C_UF4_3.CSV" "@A_SS_50C_UF5_1.CSV" "@A_SS_50C_UF5_2.CSV"
## [61] "@A_SS_50C_UF5_3.CSV" "@B_50C_LF1_1.CSV" "@B_50C_LF1_2.CSV"
## [64] "@B_50C_LF1_3.CSV" "@B_50C_LF2_1.CSV" "@B_50C_LF2_2.CSV"
## [67] "@B_50C_LF2_3.CSV" "@B_50C_LF3_1.CSV" "@B_50C_LF3_2.CSV"
## [70] "@B_50C_LF3_3.CSV" "@B_50C_LF4_1.CSV" "@B_50C_LF4_2.CSV"
## [73] "@B_50C_LF4_3.CSV" "@B_50C_LF5_1.CSV" "@B_50C_LF5_2.CSV"
## [76] "@B_50C_LF5_3.CSV" "@B_50C_UF1_1.CSV" "@B_50C_UF1_2.CSV"
## [79] "@B_50C_UF1_3.CSV" "@B_50C_UF2_1.CSV" "@B_50C_UF2_2.CSV"
## [82] "@B_50C_UF2_3.CSV" "@B_50C_UF3_1.CSV" "@B_50C_UF3_2.CSV"
## [85] "@B_50C_UF3_3.CSV" "@B_50C_UF4_1.CSV" "@B_50C_UF4_2.CSV"
## [88] "@B_50C_UF4_3.CSV" "@B_50C_UF5_1.CSV" "@B_50C_UF5_2.CSV"
## [91] "@B_50C_UF5_3.CSV" "@B_SS_50C_LF1_1.CSV" "@B_SS_50C_LF1_2.CSV"
## [94] "@B_SS_50C_LF1_3.CSV" "@B_SS_50C_LF2_1.CSV" "@B_SS_50C_LF2_2.CSV"
## [97] "@B_SS_50C_LF2_3.CSV" "@B_SS_50C_LF3_1.CSV" "@B_SS_50C_LF3_2.CSV"
## [100] "@B_SS_50C_LF3_3.CSV" "@B_SS_50C_LF4_1.CSV" "@B_SS_50C_LF4_2.CSV"
## [103] "@B_SS_50C_LF4_3.CSV" "@B_SS_50C_LF5_1.CSV" "@B_SS_50C_LF5_2.CSV"
## [106] "@B_SS_50C_LF5_3.CSV" "@B_SS_50C_UF1_1.CSV" "@B_SS_50C_UF1_2.CSV"
## [109] "@B_SS_50C_UF1_3.CSV" "@B_SS_50C_UF2_1.CSV" "@B_SS_50C_UF2_2.CSV"
## [112] "@B_SS_50C_UF2_3.CSV" "@B_SS_50C_UF3_1.CSV" "@B_SS_50C_UF3_2.CSV"
## [115] "@B_SS_50C_UF3_3.CSV" "@B_SS_50C_UF4_1.CSV" "@B_SS_50C_UF4_2.CSV"
## [118] "@B_SS_50C_UF4_3.CSV" "@B_SS_50C_UF5_1.CSV" "@B_SS_50C_UF5_2.CSV"
## [121] "@B_SS_50C_UF5_3.CSV"Revisamos si los triplicados estan completos:
(length(nombres)-1)/3## [1] 40Hasta el momento, hay un espectro repetido:
nombres[46]## [1] "@A_SS_50C_UF1_1(1).CSV"Vamos a continuar sin este espectro por el momento:
nombres <- nombres[-46]Procedemos a utilizr la función read.csv para leer los
archivos .CSV y crear objetos dentro de R que van a
contener su información.
Ya que nombres es una lista, podemos utilizar la función
de ciclo lapply para aplicar de manera iterativa la función
read.csv a todos los elementos de la lista
nombres, guardando el resultado de cada iteración en cada
elemento de la lista data.
data <- lapply(nombres, read.csv, header = FALSE)data es un objeto de clase list que
contiene en cada posición, un data.frame con la información
de los .csv. Es decir, una columna con los números de onda,
y otra con las absorbancias:
str(data[[1]])## 'data.frame': 1869 obs. of 2 variables:
## $ V1: num 399 401 403 405 407 ...
## $ V2: num 0.0129 0.0122 0.0122 0.0112 0.0106 ...En este punto, podemos guardar la información de los números de onda
en un objeto llamado wavenumber
wavenumber <- data[[1]][,1]En el chunk anterior guardamos lo que hay en la primera columna
[,1] del primer elemento de data:
data[[1]]
Ahora, podemos volver a aplicar la función lapply, en
este caso a la lista data para extraer solo la segunda
columna, que es la que nos interesa para construir nuestro data.frame
para análisis quimiométrico.
data2 <- lapply(data, '[', 2)Al usar la función lapply, de nuevo creamos una lista.
podemos convertir la lista data2 en un
data.frame mediante la función
as.data.frame.
data2 <- as.data.frame(data2)
dim(data2)## [1] 1869 120names(data2) <- nombresAhora tenemos un data.frame con 1869 filas y 120
columnas, los triplicados de 40 muestras.
Ahora podemos proceder a crear la transpuesta de este data.frame
utilizando la función t():
datat <- t(data2)datat es un objeto de clase matrix, en
cuyos rownames contiene los nombres de las muestras, como
vienen en los nombres de los archivos .CSV.
Ahora, recuperamos la información almacenada en el objeto
wavenumber para nombrar las columnas, o variables del
objeto datat
colnames(datat) <- wavenumberAhora, para encontrar los espectros a los cuales les vamos a calcular el promedio, utilizamos una matriz de búsqueda, conteniendo los niveles que diferencian las muestras, en este caso son:
- Zona de muestreo:
UF1UF2UF3UF4UF5yLF1LF2LF3LF4LF5. - tiempo de molienda:
AyB. - con o sin solución salina
SSy.
set.seed(0)
searchGrid <- expand.grid(Sampling = c('UF1',
'UF2',
'UF3',
'UF4',
'UF5',
'LF1',
'LF2',
'LF3',
'LF4',
'LF5'),
millingTime = c('A',
'B'),
brine = c('SS',
''))
# knitr::kable(searchGrid)En la matriz searchGrid tenemos un total de 40
combinaciones para los nombres de las muestras, según las distintas
variaciones de los tratamientos.
Vamos a utilizar esta matriz, para buscar entre todos los nombres de todos los espectros, y encontrar los tres triplicados para cada variación del experimento:
Primero, vamos a crear un vector que va a contener los nuevos nombres que identifican cada variación del experimento:
rownames.prom <- character(40)Luego, mediante un ciclo for podemos guardar en cada
posición de rownames.prom la unión de las tres variables
para cada experimento:
for(i in 1:40){
rownames.prom[i] <- paste(searchGrid[i,1],
searchGrid[i,2],
searchGrid[i,3])
}
rownames.prom## [1] "UF1 A SS" "UF2 A SS" "UF3 A SS" "UF4 A SS" "UF5 A SS" "LF1 A SS"
## [7] "LF2 A SS" "LF3 A SS" "LF4 A SS" "LF5 A SS" "UF1 B SS" "UF2 B SS"
## [13] "UF3 B SS" "UF4 B SS" "UF5 B SS" "LF1 B SS" "LF2 B SS" "LF3 B SS"
## [19] "LF4 B SS" "LF5 B SS" "UF1 A " "UF2 A " "UF3 A " "UF4 A "
## [25] "UF5 A " "LF1 A " "LF2 A " "LF3 A " "LF4 A " "LF5 A "
## [31] "UF1 B " "UF2 B " "UF3 B " "UF4 B " "UF5 B " "LF1 B "
## [37] "LF2 B " "LF3 B " "LF4 B " "LF5 B "Ahora, creamos la matriz en la cual vamos a guardar los resultados el promedio de cada triplicado:
prom <- matrix(,
ncol = ncol(datat),
nrow = nrow(datat)/3)Ahora, vamos a crear una lista llamada index. En esta
lista, vamos a guardar las 3 posiciones dentro de la lista
nombres que van a coincidir con las expresiones regulares
dictadas por searchGrid para cada variación del
experimento.
Por ejemplo, para UF1 A SS, la posición correspondiente
en la lista index va a contener 3 elementos que nos dicen
las posiciones para los espectros UF1 A SS 1,
UF1 A SS 2, y UF1 A SS 3 dentro de los nombres
de las filas de datat.
index <- vector('list',40)Ahora, vamos a utilizar una función de búsqueda, anidada con una función de comparación.
greplfunción de búsqueda que utiliza expresiones regulares (variables según la búsqueda) junto con caracteres (provenientes desearchGrid) como inputwhichfunción de comparación que nos dice si se cumple o no una condición, pasa como argumento el output degreplque nos dice si sí o no coincide el elemento en cuestión con el criterio de búsqueda.
for(i in 1:40){
index[[i]] <- c(
which(
grepl(
paste0('(?=.*',
as.character(searchGrid[i,1]),
')(?=.*',
as.character(searchGrid[i,2]),
')(?=.*',
as.character(searchGrid[i,3]),
')'
),
rownames(datat), perl = T )
)
)
}podemos confirmar los resultados antes de calcular. Por ejemplo, para
la primera fila de searchGrid :
searchGrid[1,]## Sampling millingTime brine
## 1 UF1 A SSYa tenemos ubicados los espectros que corresponden a esta combinación:
index[[1]]## [1] 46 47 48confirmando en los nombres de las muestras de datat:
rownames(datat)[index[[1]]]## [1] "@A_SS_50C_UF1_1.CSV" "@A_SS_50C_UF1_2.CSV" "@A_SS_50C_UF1_3.CSV"Corresponde a los tres triplicados de la muestra UF1, la
cual tuvo un tiempo de molienda de 2 a 4 minutos, y una extracción con
solución salina.
Podemos proceder a calcular los promedios:
for(j in 1:length(colnames(datat))){
for(i in 1:40){
prom[i,j] <- mean(c(datat[index[[i]][1],j],
datat[index[[i]][2],j],
datat[index[[i]][3],j]
) )
}
}
rownames(prom) <- rownames.prom
colnames(prom) <- wavenumberUna vez calculados, podemos graficar los promedios:
Podemos diferenciar los espectros de las muestras a las cuales se les realizo la extracción con solución salina mediante una revisión
NaCl <- vector('character', length(rownames(prom)))
for(i in 1:40){
if(grepl('SS', rownames(prom)[i],perl =T)){NaCl[i]<- 'red'}
else{NaCl[i]<- 'black'}
}
indexSS <- vector('logical',40)
indexSS<-which(grepl('SS', rownames(prom),perl =T))for(i in 1:40){
plot(wavenumber,
prom[i,],
xlab = 'wave number cm-1',
ylab = 'absorbance a.u.',
xlim = c(3900,400),
ylim = c(0,0.12),
type = 'l',
col = NaCl[i])
par(new = T)
}
Podemos seleccionar este rango , de 1700 a 400 \(cm^{-1}\) creando un nuevo objeto para los promedios y para los números de ondaÑ
colnames(prom)[c(1,676)]## [1] "399.1989" "1700.935"Rango1 <- prom[,1:676]
wn.Rango1 <- wavenumber[1:676]Rango2 <- Rango1[, -c(1:163)]
wn.Rango2 <- wn.Rango1[-c(1:163)] # Quitando desde 790 hacia abajoPara eliminar la variabilidad por ruido instrumental, podemos
calcular y restar la línea base para cada uno de los promedios,
utilizando el método rubberband disponible en el paquete de
R hyperSpec, desarrollado para quimiometría aplicada a
imágenes hiper-espectrales.
library(hyperSpec)
spc <- new('hyperSpec', spc = Rango2, wavelength = wn.Rango2)
bend <- 0.1 * wl.eval (spc, function (x) x^6+x^5+x^4+x^3+x^2, normalize.wl=normalize01)
bl <- spc.rubberband (spc+bend, noise = 1e-4, df=20)-bend
suma <- spc+bend
spc3 <- spc - bl
plot(spc, wl.reverse = TRUE)
plot(bl, add=TRUE, col=2,wl.reverse = TRUE)
plot(suma,wl.reverse = TRUE)
plot(bend, add=TRUE, col=2,wl.reverse = TRUE)
plot(spc3,wl.reverse = TRUE)
corregido <- as.data.frame(spc3[1:40])
corregido2 <- as.data.frame(corregido[,1])Hasta este punto se conserva una parte de la variabilidad por ruido instrumental en el rango de 790 a 400 cm-1. Dado a que en el rango de 1700 a 790 cm-1 tenemos información de enlaces siloxano en dos bandas, podemos utilizar este rango para el analisis multivariado:
Podemos graficar los espectros corregidos para este experimento:
for(i in 1:length(rownames(corregido2))){
plot(wn.Rango2,
corregido2[i,],
xlab = 'wave number cm-1',
ylab = 'absorbance a.u.',
xlim = c(1700,790),
ylim = c(0,0.06),
type = 'l',
col = NaCl[i])
par(new = T)
}
legend('topleft',
c('50C Colombia SS','50C Colombia'),
col = c('red','black'),
lty = 1)
pca.bl <- prcomp(corregido2, scale =T, center = T)
vp.bl <- (pca.bl$sdev)^2
varianza.bl <- round(vp.bl/sum(vp.bl)*100, 2)
varianza.bl## [1] 63.88 8.64 8.28 6.40 3.20 1.89 1.55 1.25 0.76 0.57 0.49 0.43
## [13] 0.33 0.30 0.22 0.20 0.18 0.16 0.14 0.12 0.11 0.10 0.09 0.08
## [25] 0.08 0.07 0.06 0.06 0.05 0.05 0.04 0.04 0.04 0.03 0.03 0.02
## [37] 0.02 0.02 0.01 0.00coord.bl <- pca.bl$x
loadings <- pca.bl$rotationlibrary(RColorBrewer)
for(i in 1:2){
plot(as.numeric(colnames(corregido2)),
loadings[,i],
type = 'l',
ylim = c(-0.15,0.2),
xlim = c(1700,790),
ylab = 'PC Loadings',
xlab = expression(paste("Wave number (cm"^"-1",")")),
col= brewer.pal(3,'Dark2')[i],
)
par(new =T)
}
abline(h=0)
legend('topleft',
c('PC1','PC2'),
col = brewer.pal(3,'Dark2')[1:2],
lty = 1,
lw = 2)
Derivatización de espectros
library(prospectr)## [34mprospectr version 0.2.5 -- [39m'antilla'## [34mcheck the github repository at: https://github.com/l-ramirez-lopez/prospectr/[39mw2 <- 40
sg <- matrix(ncol = ncol(corregido2)-w2, nrow= nrow(corregido2))
sg <- as.data.frame(sg)
# for (i in 1:31){
sg<- savitzkyGolay(X = corregido2
,m = 2,
p = 3,
w = w2+1,
delta.wav=2)
# }
# colnames(sg) <- colnames( savitzkyGolay(X = corregido2[1,]
# ,m = 2,
# p = 3,
# w = w2+1,
# delta.wav=2))
rownames(sg) <- rownames(corregido2)
# win.graph()
for(i in 1:length(rownames(sg))){
plot(as.numeric(colnames(sg)),
sg[i,],
xlab = 'wave number cm-1',
ylab = '2nd derivative of abs',
xlim = c(1700,790),
ylim = c(-5e-05,5e-05),
type = 'l',
col = NaCl[i])
par(new = T)
}
legend('topleft',
c('50C Colombia SS','50C Colombia'),
col = c('red','black'),
lty = 1)
Análisis exploratorio
Cambio de nombres para hca
nombres.hca <- vector('character', 40)
for (i in 1:40){
if(grepl('(?=.*UF)(?=.*A)((?!SS).)*$',rownames(sg)[i],perl=T)){nombres.hca[i]<- '@50C-Colombia 2-4 min'}else{
if(grepl('(?=.*UF)(?=.*B)((?!SS).)*$',rownames(sg)[i],perl=T)){nombres.hca[i]<- '@50C-Colombia 6-8 min'}else{
if(grepl('(?=.*UF)(?=.*A)(?=.*SS)',rownames(sg)[i],perl=T)){nombres.hca[i]<- '@50C-Colombia 2-4 min SS'}else{
if(grepl('(?=.*UF)(?=.*B)(?=.*SS)',rownames(sg)[i],perl=T)){nombres.hca[i]<- '@50C-Colombia 6-8 min SS'}else{
if(grepl('(?=.*LF)(?=.*A)((?!SS).)*$',rownames(sg)[i],perl=T)){nombres.hca[i]<- '50C-Colombia 2-4 min'}else{
if(grepl('(?=.*LF)(?=.*B)((?!SS).)*$',rownames(sg)[i],perl=T)){nombres.hca[i]<- '50C-Colombia 6-8 min'}else{
if(grepl('(?=.*LF)(?=.*A)(?=.*SS)',rownames(sg)[i],perl=T)){nombres.hca[i]<- '50C-Colombia 2-4 min SS'}else{
if(grepl('(?=.*LF)(?=.*B)(?=.*SS)',rownames(sg)[i],perl=T)){nombres.hca[i]<- '50C-Colombia 6-8 min SS'}else{
}
}
}
}
}
}
}
}
}compare <- data.frame(antes = rownames(sg), ahora = nombres.hca)rownames(sg) <- nombres.hcafraction <- vector('integer',40)
for(i in 1:40){
if(grepl('(?=.*@)',rownames(sg)[i],perl=T)){fraction[i] <- 1}else{fraction[i] <- 2}
}
# Upper 1
# Lower 2HCA
df <- sg
library(factoextra)## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBares.hk <-hkmeans(df,
5,
hc.metric = 'minkowski')
fviz_dend(res.hk, cex = 0.6, palette = "jco",
rect = TRUE, rect_border = "jco", rect_fill = TRUE, ylim = c(-5e-05,2e-04))## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
pca.bl.sg <- prcomp(sg, scale =T, center = T)
vp.bl.sg <- (pca.bl.sg$sdev)^2
varianza.bl.sg <- round(vp.bl.sg/sum(vp.bl.sg)*100, 2)
varianza.bl.sg## [1] 70.16 14.36 7.84 3.11 1.25 0.97 0.88 0.46 0.28 0.18 0.13 0.10
## [13] 0.06 0.05 0.03 0.03 0.03 0.02 0.01 0.01 0.01 0.01 0.00 0.00
## [25] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## [37] 0.00 0.00 0.00 0.00coord.bl.sg <- pca.bl.sg$xplot(coord.bl.sg[,1],
coord.bl.sg[,2],
pch = 19,
xlab = paste('PC 1 -',as.character(varianza.bl.sg[1]), '%'),
ylab = paste('PC 2 -',as.character(varianza.bl.sg[2]), '%'),
col = res.hk$cluster,
xlim = c(-20,90)
)
abline(h = 0, v= 0, lty = 2)
x<- coord.bl.sg[,1]
y<- coord.bl.sg[,2]
z<- coord.bl.sg[,3]
x1<- x[res.hk$cluster == 1]
y1<- y[res.hk$cluster == 1]
z1<- z[res.hk$cluster == 1]
library(car)## Loading required package: carDataellipse(c(mean(x1),mean(y1)),
cov(cbind(x1,y1)),
radius= sqrt(qchisq(.95, df =2)),
col=1,
center.pch=FALSE,
add=TRUE,
# levels=0.95,
fill=TRUE,
lwd=0.8,
)
x2<- x[res.hk$cluster == 2]
y2<- y[res.hk$cluster == 2]
z2<- z[res.hk$cluster == 2]
ellipse(c(mean(x2),mean(y2)),
cov(cbind(x2,y2)),
radius= sqrt(qchisq(.95, df =2)),
col=2,
center.pch=FALSE,
add=TRUE,
# levels=0.95,
fill=TRUE,
lwd=0.8,
)
x3<- x[res.hk$cluster == 3]
y3<- y[res.hk$cluster == 3]
z3<- z[res.hk$cluster == 3]
ellipse(c(mean(x3),mean(y3)),
cov(cbind(x3,y3)),
radius= sqrt(qchisq(.95, df =2)),
col=3,
center.pch=FALSE,
add=TRUE,
# levels=0.95,
fill=TRUE,
lwd=0.8,
) 
# text(coord.bl.sg[,1],
# coord.bl.sg[,2],
# rownames(coord.bl.sg),
# pos = 4,
# cex = 0.5)x <- coord.bl.sg[,1]
y <- coord.bl.sg[,2]
z <- coord.bl.sg[,3]
x1 <- x[NaCl=='red']
y1 <- y[NaCl=='red']
z1 <- z[NaCl=='red']
x2 <- x[NaCl=='black']
y2 <- y[NaCl=='black']
z2 <- z[NaCl=='black']
plot(coord.bl.sg[,1],
coord.bl.sg[,2],
pch = 19,
xlab = paste('PC 1 -',as.character(varianza.bl.sg[1]), '%'),
ylab = paste('PC 2 -',as.character(varianza.bl.sg[2]), '%'),
col = NaCl
)
abline(h = 0, v= 0, lty = 2)
library(car)
ellipse(c(mean(x1),mean(y1)),
cov(cbind(x1,y1)),
radius= sqrt(qchisq(.95, df =2)),
col="red",
center.pch=FALSE,
add=TRUE,
# levels=0.95,
fill=TRUE,
lwd=0.8,
) 
# ellipse(c(mean(x2),mean(y2)),
# cov(cbind(x2,y2)),
# radius= sqrt(qnorm(.95)),
# col="black",
# center.pch=FALSE,
# add=TRUE,
# # levels=0.95,
# fill=TRUE,
# lwd=0.8,
#
# ) library(mdatools)##
## Attaching package: 'mdatools'## The following object is masked from 'package:car':
##
## ellipsemdaPCA <-pca(corregido2, 5, center = TRUE, scale = TRUE)
summary(mdaPCA)##
## Summary for PCA model (class pca)
## Type of limits: ddmoments
## Alpha: 0.05
## Gamma: 0.01
##
## Eigenvals Expvar Cumexpvar Nq Nh
## Comp 1 327.705 63.88 63.88 1 1
## Comp 2 44.304 8.64 72.52 2 1
## Comp 3 42.459 8.28 80.79 4 1
## Comp 4 32.839 6.40 87.19 6 1
## Comp 5 16.402 3.20 90.39 5 1plot(mdaPCA, show.labels = T)
# win.graph()
plotScores(mdaPCA$res$cal, show.labels = TRUE)
Ya que mediante PCA para estas muestras en específico, no se obserban
agrupamientos notorios, podemos exportar estos resultados y compararlos
con los resultados de los experimentos anteriores:
write.table(prom, 'SSexp.txt' , sep = '\t')