Clase 6 Multivariado covariaciones, supuestos de normalidad y MANCOVA
Javier Borbon
2024-03-21
#en el modelado tengo tres respuestas, en 10 trt, cual incentiva mas el crc de las bacterias?
#Hacer el analisis hipotesis x H2 x H3 X H4, suele estar mal hacer el analisis variable x variable
#Cargar datos de papa
library(readxl)
Datos_papa <- read_excel("C:/2024/Semestre III/Metodos Multivariados/Datos_papa.xlsx")
df=Datos_papa#vamos a hacer una matrix de correlaciones, para entender xk el analisis debe ser covariado
library(corrplot)## corrplot 0.92 loadedcl = cor(df[ ,1:3]) #una coma con espacio vacio antes, en el data.frame tengo filas y columnas, si quiero leer todas las filas, dejo el espacio vacio, pero lease las columnas 1 a 3, la 4 son los trt
corrplot(cl, method = "n",
number.cex = 3)
library(ggplot2)
ggplot(df, aes(x = LR, y = VR))+
geom_point()+
geom_boxplot(aes(fill = trt))
ggplot(df, aes(x = LR, y = PSR))+
geom_point()+
geom_boxplot(aes(fill = trt))
ggplot(df, aes(x = VR, y = PSR))+
geom_point()+
geom_boxplot(aes(fill = trt))
#relacionamos volumen y largo de raiz, y se ve una tendencia entre los puntos, es posible que tenga correlaciones
#el control se me puede dibujar en otra zona, mostrando que este se porto mas malito
#en el grafico de PSR x VR el largo de raices puede estar en el ancho de las cajas#hagamos la nube de puntos
library(rgl)
cols=c("brown","red","blue","coral1", "green", "black", "cyan", "yellow","deeppink", "orange")
col = cols[as.factor(df$trt)]
legend3d("topright",
legend = unique(df$trt), pch=19,
col=cols[as.factor(unique(df$trt))],
cex=0.9, inset=c(0.03))
plot3d(df$LR ,df$VR ,df$PSR, type="s",
col=col, radius = .1, data=df,
size = 10)#descripcion de la estadistica de los datos
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(psych)##
## Attaching package: 'psych'## The following objects are masked from 'package:ggplot2':
##
## %+%, alphapsych::describeBy(x = df[ , 1:3],
group = df$trt)##
## Descriptive statistics by group
## group: BC
## vars n mean sd median trimmed mad min max range skew kurtosis
## LR 1 12 29.27 1.83 29.67 29.31 1.27 26.03 32.09 6.07 -0.46 -0.97
## VR 2 12 5.89 0.94 6.02 5.90 0.94 4.32 7.34 3.02 -0.29 -1.20
## PSR 3 12 5.23 1.02 5.34 5.29 1.14 3.33 6.47 3.13 -0.53 -1.12
## se
## LR 0.53
## VR 0.27
## PSR 0.29
## ------------------------------------------------------------
## group: BM
## vars n mean sd median trimmed mad min max range skew kurtosis
## LR 1 12 31.49 2.04 32.07 31.71 1.67 26.69 34.07 7.39 -0.89 0.01
## VR 2 12 6.00 0.75 6.00 6.04 0.53 4.48 7.20 2.71 -0.25 -0.65
## PSR 3 12 5.52 0.89 5.75 5.51 0.29 3.88 7.19 3.31 -0.35 -0.32
## se
## LR 0.59
## VR 0.22
## PSR 0.26
## ------------------------------------------------------------
## group: BMM
## vars n mean sd median trimmed mad min max range skew kurtosis se
## LR 1 12 31.18 1.26 31.04 31.17 0.85 28.64 33.74 5.10 0.09 -0.07 0.36
## VR 2 12 5.88 0.56 5.79 5.88 0.47 4.87 6.93 2.06 0.24 -0.61 0.16
## PSR 3 12 5.04 0.91 4.96 4.94 0.80 3.85 7.22 3.37 0.89 0.21 0.26
## ------------------------------------------------------------
## group: EA
## vars n mean sd median trimmed mad min max range skew kurtosis
## LR 1 12 32.57 1.49 32.80 32.52 1.27 30.31 35.33 5.02 0.07 -0.94
## VR 2 12 6.29 1.07 6.39 6.25 0.91 4.49 8.43 3.94 0.08 -0.57
## PSR 3 12 5.53 1.02 5.87 5.56 1.11 3.78 6.99 3.21 -0.34 -1.29
## se
## LR 0.43
## VR 0.31
## PSR 0.29
## ------------------------------------------------------------
## group: M
## vars n mean sd median trimmed mad min max range skew kurtosis
## LR 1 12 32.74 0.96 32.88 32.77 1.00 31.06 34.11 3.06 -0.27 -1.30
## VR 2 12 6.48 0.57 6.28 6.44 0.54 5.68 7.62 1.94 0.47 -1.04
## PSR 3 12 6.42 0.82 6.44 6.36 0.68 5.06 8.36 3.30 0.70 0.34
## se
## LR 0.28
## VR 0.17
## PSR 0.24
## ------------------------------------------------------------
## group: MT
## vars n mean sd median trimmed mad min max range skew kurtosis se
## LR 1 12 31.79 1.15 31.91 31.73 1.43 30.37 33.78 3.41 0.11 -1.52 0.33
## VR 2 12 6.07 0.58 6.13 6.05 0.62 5.22 7.09 1.88 0.07 -1.28 0.17
## PSR 3 12 5.78 0.79 5.73 5.76 0.37 4.34 7.39 3.06 0.22 -0.30 0.23
## ------------------------------------------------------------
## group: OS
## vars n mean sd median trimmed mad min max range skew kurtosis
## LR 1 12 32.47 1.30 32.74 32.62 1.41 29.65 33.81 4.16 -0.73 -0.65
## VR 2 12 6.22 0.67 6.32 6.19 0.62 5.31 7.41 2.11 0.00 -1.29
## PSR 3 12 5.48 0.80 5.44 5.43 0.56 4.26 7.18 2.92 0.47 -0.49
## se
## LR 0.37
## VR 0.19
## PSR 0.23
## ------------------------------------------------------------
## group: PF
## vars n mean sd median trimmed mad min max range skew kurtosis se
## LR 1 12 29.16 1.23 29.45 29.11 0.41 26.82 31.97 5.14 0.33 0.38 0.36
## VR 2 12 6.00 0.67 5.91 5.94 0.57 5.13 7.44 2.30 0.63 -0.57 0.19
## PSR 3 12 5.23 0.74 5.14 5.23 0.88 3.93 6.56 2.62 0.01 -1.02 0.21
## ------------------------------------------------------------
## group: PFM
## vars n mean sd median trimmed mad min max range skew kurtosis
## LR 1 12 28.86 1.56 28.26 28.73 0.85 26.77 32.18 5.41 0.69 -0.73
## VR 2 12 6.07 0.54 6.17 6.10 0.60 5.07 6.79 1.73 -0.40 -1.21
## PSR 3 12 5.08 0.65 4.93 5.08 0.74 4.14 5.98 1.83 0.13 -1.58
## se
## LR 0.45
## VR 0.16
## PSR 0.19
## ------------------------------------------------------------
## group: T0
## vars n mean sd median trimmed mad min max range skew kurtosis
## LR 1 12 22.13 0.92 22.11 22.12 1.24 20.85 23.56 2.71 0.15 -1.59
## VR 2 12 5.07 0.61 5.11 5.07 0.70 4.06 6.12 2.06 0.02 -1.19
## PSR 3 12 4.90 0.79 5.17 4.95 0.50 3.38 5.90 2.53 -0.73 -0.89
## se
## LR 0.26
## VR 0.18
## PSR 0.23#esta la estadistica por cada trt
#Kurtosis,
#medidas de dispersión, desviacion estadar, media, rango
#localizacion: min, max, cuartiles, percentiles
#forma: forma de la curva: kurtosis y simetria,
#skew, simetria, hasta que punto tengo datos simetricos, que es lo normal
##meso curtica, planicurtica, platicurtica #libreria de datos
library(PerformanceAnalytics)## Loading required package: xts## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric##
## ################################### WARNING ###################################
## # We noticed you have dplyr installed. The dplyr lag() function breaks how #
## # base R's lag() function is supposed to work, which breaks lag(my_xts). #
## # #
## # Calls to lag(my_xts) that you enter or source() into this session won't #
## # work correctly. #
## # #
## # All package code is unaffected because it is protected by the R namespace #
## # mechanism. #
## # #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning. #
## # #
## # You can use stats::lag() to make sure you're not using dplyr::lag(), or you #
## # can add conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop #
## # dplyr from breaking base R's lag() function. #
## ################################### WARNING #####################################
## Attaching package: 'xts'## The following objects are masked from 'package:dplyr':
##
## first, last##
## Attaching package: 'PerformanceAnalytics'## The following object is masked from 'package:graphics':
##
## legendPerformanceAnalytics::chart.Correlation(R = df[ ,1:3])## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
##grafica 1 pareciera ser bimodal, la grafica de la mitad se ve normal
#el grafico me dice que es evidente que las variables estan relacionadas
#aunque se observa la normalidad de VR, debo hacer analisis multivariado
#PERO ESTAS GRAFICAS ME ESTAN AGRUPANDO DATOS DE LOS 10 TRT, aca veo el comportamiento de todos los datos
#Si hubiese hecho los graficos por trt, tengo un problema, solo 12 resultados por trt, tengo poca exactitud
##trt por color
plot(df$LR, df$VR, pch = 16,
col = as.factor(df$trt))
#veo colores en sitios diferentes, se diferencia el control a la izq, una agrupacion entre negro y gris, y los demás a la derecha
#tendría 3 trt, o agrupación
#hagamoslos en las otras variables
plot(df$PSR, df$VR, pch = 16,
col = as.factor(df$trt))
#aca no se ve separacion de control
plot(df$LR, df$PSR, pch = 16,
col = as.factor(df$trt))
#aca vuelven a verse 3 trt (3 xk los agrupe)#hagamos un pipe
#partirme por trt
df |>
group_by(trt) |>
summarise(Media_LR = mean(LR),
Desv_LR = sd(LR),
CV_LR = Desv_LR/Media_LR,
Media_VR = mean(VR),
Desv_VR = sd(VR),
CV_VR = Desv_VR/Media_VR,
Media_PSR = mean(PSR),
Desv_PSR = sd(PSR),
CV_PSR = Desv_PSR/Media_PSR) ## # A tibble: 10 × 10
## trt Media_LR Desv_LR CV_LR Media_VR Desv_VR CV_VR Media_PSR Desv_PSR
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 BC 29.3 1.83 0.0624 5.89 0.945 0.160 5.23 1.02
## 2 BM 31.5 2.04 0.0646 6.00 0.748 0.125 5.52 0.891
## 3 BMM 31.2 1.26 0.0405 5.88 0.560 0.0953 5.04 0.908
## 4 EA 32.6 1.49 0.0459 6.29 1.07 0.170 5.53 1.02
## 5 M 32.7 0.964 0.0295 6.48 0.573 0.0885 6.42 0.822
## 6 MT 31.8 1.15 0.0362 6.07 0.577 0.0952 5.78 0.795
## 7 OS 32.5 1.30 0.0400 6.22 0.666 0.107 5.48 0.805
## 8 PF 29.2 1.23 0.0423 6.00 0.670 0.112 5.23 0.744
## 9 PFM 28.9 1.56 0.0541 6.07 0.537 0.0884 5.08 0.651
## 10 T0 22.1 0.916 0.0414 5.07 0.607 0.120 4.90 0.788
## # ℹ 1 more variable: CV_PSR <dbl>#hagamos la parte descriptiva
#PSR, en micorrizas se lleva los 3 promedios mas altos, pareciera que micorrizas lo induce
#Medir, tabular y
#1 compenente descriptivo, 2 componente inferencial (de analisis)Modelo \[Y_{ij} = \mu + \tau_i + \epsilon_{ij}\] \[i= 1,...,10\] \[j= 1,...,12\]
Analisis de varianza
#Diseño experimental: Factoral simple en arreglo completamente al azar (no tengo bloques)
#modelo si hago el analisis variable por variable - barra inversa\ -
#Modelo Largo raices
mod_LR = aov(formula = LR ~ trt, data = df)
summary(mod_LR)## Df Sum Sq Mean Sq F value Pr(>F)
## trt 9 1093.8 121.5 60.63 <2e-16 ***
## Residuals 110 220.5 2.0
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#el largo de las raices, es funcion de los tratamientos, esto se hace con virgulilla
#Analisis de varianza alt 126#en la tabla anterior, hubo alta significancia, P= 2x10 a la -16, esto se aproxima a 0 xk es el numero mas bajito de R, deberiamos ponerlo similar o igual a 0
#rechazo H0, pero error, tengo 3 hipotesis, entonces el valor P ya no es 0.05
#este ANOVA me dice que rechazo, entonces al menos 1 de los valores es diferente
#Ahora cambio de variable, lo llamo modelo 2 (VR) y lo corro
mod_VR = aov(formula = VR ~ trt, data = df)
summary(mod_VR)## Df Sum Sq Mean Sq F value Pr(>F)
## trt 9 15.11 1.6786 3.277 0.00143 **
## Residuals 110 56.34 0.5122
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod_PSR = aov(formula = PSR ~ trt, data = df)
summary(mod_PSR)## Df Sum Sq Mean Sq F value Pr(>F)
## trt 9 21.12 2.346 3.241 0.00158 **
## Residuals 110 79.64 0.724
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#En los 3 modelos parece que los trt difieren, pero no sabemos quien es el diferente
#estoy omitiendo las correlaciones
#estoy omitiendo el p valor ajustado por la relacion entre variables#supuestos de normalidad
#metodo parametrico que poco se enseña: kruskal wallis
#prueba shapiro-wilks
shapiro.test(mod_LR$residuals)##
## Shapiro-Wilk normality test
##
## data: mod_LR$residuals
## W = 0.98599, p-value = 0.2515shapiro.test(mod_VR$residuals)##
## Shapiro-Wilk normality test
##
## data: mod_VR$residuals
## W = 0.99371, p-value = 0.87shapiro.test(mod_PSR$residuals)##
## Shapiro-Wilk normality test
##
## data: mod_PSR$residuals
## W = 0.99061, p-value = 0.5896#la diferencia entre lo ajustado (proyectado vs real?) es el residual
#las tres variables ya casi son normales... PEEEERO la normalidad variable por variable no garantizar que sea normal, pero tengo que evaluarlo trivariantemente
#si hay normalidad variable x variable, no significa que la normalidad bivariante o trivariante tmbn lo sea
#bivariante, la grafica parece una espinilla, trivariante ya ni se puede graficar#otro supuesto, homostacilidad, que los trt tengan la misma variabilidad
#hagamos un boxplot
#prueba de Homocedasticidad, igualdad de varianza
boxplot(LR ~ trt, df)
#si las varianzas no son iguales?
#comparo las medias, pero la varianza de uno es 5 veces mas grande que la otra,
#me significa que las varianzas son desiguales, las medias no son comparables
#si varia mucho, ni siquiera puedo creer en su promedio, ni los puedo comparar
#si quiero comparar, debo iniciar con condiciones inicales
#estas grandes varianzas, las estandarizo, penalizo, se hacen promedios ponderados, se ajusta para que haya razonabilidad en las comparaciones
#dos dupuestos que se volvieron costumbre: normalidad y analisis de varianza
#tmbn el efecto borde me daña mi estadsitica: poligonos de Voronoi o teselacion: se generan datos virtuales para completar 1 a 1 dato y poder graficar o comparar, pero hay mucho error aqui#test de barlet para verificar homocedasticidad
bartlett.test(df$LR, df$trt)##
## Bartlett test of homogeneity of variances
##
## data: df$LR and df$trt
## Bartlett's K-squared = 12.677, df = 9, p-value = 0.1778bartlett.test(df$VR, df$trt)##
## Bartlett test of homogeneity of variances
##
## data: df$VR and df$trt
## Bartlett's K-squared = 11.685, df = 9, p-value = 0.2316bartlett.test(df$PSR, df$trt)##
## Bartlett test of homogeneity of variances
##
## data: df$PSR and df$trt
## Bartlett's K-squared = 3.6044, df = 9, p-value = 0.9355#errores de la estadistica, hacer analisis por separado
#varianzas VR 0.2 varianzas iguales, PSR 0.9 o 90%
#entonces alguien normalmente diria que hay poca o alta homocedasticidad haría lo siguiente para encontrar que trt difiere#libreria mulcom, es para deterctar la diferencia
####el siguiente codigo no se pudo correr
#par(mai=c(1,1,1.25,1))
#library(multcomp)
#df.aov1 = aov (LR ~ as.factor(trt), df)
#tukey1 = glht(df.aov1, linfct = mcp(trt="Tukey"))
#cld(tukey1)
#plot(cld(tukey1))####se retoma desde aca
s = TukeyHSD(mod_LR)
plot(s, las = 2)
#se tiene una linea continua que es mi linea 0, los intervalos que esten en esta linea, no difieren entre silibrary(agricolae)##
## Attaching package: 'agricolae'## The following objects are masked from 'package:PerformanceAnalytics':
##
## kurtosis, skewnessout_LR = duncan.test(mod_LR, "trt")
plot(out_LR, variation = "IQR")
out_VR = duncan.test(mod_VR, "trt")
plot(out_VR, variation = "IQR")
out_PSR = duncan.test(mod_PSR, "trt")
plot(out_PSR, variation = "IQR")
#Segun estas graficas, hay varios grupos
#Aca finaliza en enfoque univariado, el que solemos hacer por tradicion, el que no deberiamos hacer#Multivariado
#analisis de varianza multivariante:MANOVA
#sitioweb para codigo en R: chatgpt y quick-R
Y = cbind(df$LR, df$VR, df$PSR) #meter en una variable todo, en este caso, las 3 variables
modelo = manova(Y ~ df$trt)
summary(modelo, test="Wilks")## Df Wilks approx F num Df den Df Pr(>F)
## df$trt 9 0.086679 15.338 27 316.06 < 2.2e-16 ***
## Residuals 110
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#aca el valor p me rechaza la igualdada de trt en todas las variables, dado p = 2x10 a la -16
#lo que se es que los trt difieren para las 3
#si tuviera mas variables, puedo poner aqui modelo = manova(Y ~ df$trt) + factor + factor para que me tome los cultivos o genotipos
#si aparece un error de singular, me esta diciendo que tengo combiancion lineal de mis variables, por ejemplo:
#metí arena, arcilla y limo, si la suma es 100, si tengo dos, la tercera sale por diferencia
##una variable es la otra multiplicado por una constante
#si el determinante me da 0, no se puede invertir...
#singularidad, tenemos variables que sobran o no estan aportando#PARENTESIS MOSTRAR SINGULARIDAD (ESTO NO SE PUEDE CORRER, SON EJEMPLOS DE COMO R DETECTA SINGULARIDAD Y NO FUNCIONA, LO PONGO EN ###)
#ml = matrix(c( 1, -1, 3, 2, -2, 6, 7, 5, 1),
# nrow = 3, ncol = 3, byrow = F )
#det(ml)
#solve(ml)
#la segunda columna es la primera multiplicada x 2#PARENTESIS MOSTRAR SINGULARIDAD
#ml = matrix(c( 1, -1, 3, 8, 4, 4, 7, 5, 1),
# nrow = 3, ncol = 3, byrow = F )
#det(ml)
#solve(ml)
#la segunda columna es la primera multiplicada x 2 #ejemplo arena arcilla limo
#tiene singularidad, no se puede resolver el problema#hagamos prueba de normalidad multivariante
library(mvnTest)## Loading required package: mvtnormmvnTest::R.test(df[ ,1:3]) ## Royston test for Multivariate Normality
##
## data : df[, 1:3]
##
## R : 32.69715
## p-value : 3.361028e-07
##
## Result : Data are not multivariate normal (sig.level = 0.05)#este p valor me dice que no hay normalidad en mis variables, y aca a llorar al monte? tenemos prueba\[H_0: \Sigma_0 = \Sigma_1 = \dots = \Sigma...{10}\]
#homogenidad de covarianzas
library(biotools)## Loading required package: MASS##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## select## ---
## biotools version 4.2boxM(df[ , 1:3], df$trt)##
## Box's M-test for Homogeneity of Covariance Matrices
##
## data: df[, 1:3]
## Chi-Sq (approx.) = 46.713, df = 54, p-value = 0.7487#cumpli las varianzas (P = 0.74), pero no la normalidad
#no deberia ponerme a mirar los datos para ver como arreglo mi normalidad#dejé el tiempo de germinación por fuera, estoy grave con la normalidad, ajustemos los datos a ver si con la covariable me sirva, probemos
df$t_germ = round(sort(rnorm(n = 120, mean = 20, sd =5)))
#tiempo de germinación#ANALISIS DE COVARIANZA UNIVARIANTE/MULTIVARIANTE
#correr un modelo involucrando la germinacion
mod_co_LR = aov(df$LR ~ df$t_germ + df$trt)
summary(mod_co_LR)## Df Sum Sq Mean Sq F value Pr(>F)
## df$t_germ 1 788.6 788.6 402.75 <2e-16 ***
## df$trt 9 312.3 34.7 17.72 <2e-16 ***
## Residuals 109 213.4 2.0
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#no voy a correr tres ancovas, #MANCOVA, aca ya metemos germinacion en el modelo
mod_co_LR2 = cbind(df$LR, df$VR, df$PSR) #meter en una variable todo, en este caso, las 3 variables
modelo2 = manova(Y ~ df$t_germ + df$trt)
summary(modelo2, test="Wilks")## Df Wilks approx F num Df den Df Pr(>F)
## df$t_germ 1 0.14923 203.343 3 107.00 < 2.2e-16 ***
## df$trt 9 0.23619 7.412 27 313.14 < 2.2e-16 ***
## Residuals 109
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#acostumbremonos primero a tener las covariables y luego los factores#probemos la normalidad
#se le prueba a las variables o a los residuales? es a los RESIDUALES
#Extraemos los residuales y le hacemos la prueba de normalidad
#xk los residuales... complementar con el audio
mod_co_LR = cbind(df$LR,df$VR,df$PSR)
modelo <- manova(Y ~ df$t_germ + df$trt )
summary(modelo, test="Wilks")## Df Wilks approx F num Df den Df Pr(>F)
## df$t_germ 1 0.14923 203.343 3 107.00 < 2.2e-16 ***
## df$trt 9 0.23619 7.412 27 313.14 < 2.2e-16 ***
## Residuals 109
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#multivariable
mvnTest::R.test(modelo$residuals)## Royston test for Multivariate Normality
##
## data : modelo$residuals
##
## R : 0.7377692
## p-value : 0.8507392
##
## Result : Data are multivariate normal (sig.level = 0.05)#p valor = 0.95, los residuales son normales