Universidad Nacional de Colombia

TALLER ANÁLISIS DE DATOS EN R

1. Verifique que tenga la versión de RStudio más reciente (1.4.1717). Si no, descárguela del link: https://www.rstudio.com/products/rstudio/download/#download

DCA

2. Organice los datos que desea analizar en R del Excel o formato en donde registró sus resultados. Para este taller usaremos primero los datos de Tasa Respiratoria del Excel adjunto. Organice los datos en un nuevo Excel de tal forma que los datos queden en columnas de tratamientos, bloques (si los hay), y repeticiones. Guarde el archivo en formato Excel “.xlsx”

Nota: Los decimales deben estar en formato de punto, NO de coma, si no el software no los podrá leer. Para cambiar esta característica en Excel, antes de pegar sus datos diríjase a Archivo>Opciones>Avanzadas>Separador decimal, y ponga “.”> Aceptar.

3. Con sus datos organizados inicie Rstudio. Abra un nuevo script: File> New File> Rscript.

4. Llame las librerías a utilizar. En este caso, para este taller, usaremos las librerías:

library(readxl) 
library(agricolae)
library(car) 
library(carData) 
library(MASS)
library(RColorBrewer) #Opcional

Nota: Si no las tiene instaladas diríjase a Tool> Install Packages, inserte los nombres de las librerías a descargar separadas por coma o espacio. El programa tardará unos minutos en instalarlos.

5. Proceda a subir los datos del Excel, dándole un nombre

TasaR<- read_excel("C:/Users/Juliana Miranda/Documents/Universidad/TESIS/RESULTADOS/R/Respiracion.xlsx")

6. Visualice un resumen de sus datos:

summary(TasaR)

      DDC             TR       
 Min.   : 3.0   Min.   :11.29  
 1st Qu.: 5.0   1st Qu.:37.16  
 Median : 7.5   Median :45.90  
 Mean   : 7.5   Mean   :45.09  
 3rd Qu.:10.0   3rd Qu.:55.06  
 Max.   :12.0   Max.   :81.45  

7. Establezca los días después de cosecha (DDC) como factor

TasaR$DDC=as.factor(TasaR$DDC)

8. Establezca el modelo para un DCA simple y genere el análisis de varianza

TRanov=lm(TasaR$TR~TasaR$DDC)
anova(TRanov)

Analysis of Variance Table

Response: TasaR$TR
          Df  Sum Sq Mean Sq F value    Pr(>F)    
TasaR$DDC  9 15042.9 1671.43  17.549 < 2.2e-16 ***
Residuals 90  8572.1   95.25                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

¿A que hace referencia cada uno de los items de la tabla? ¿qué se puede inferir de esta información? ¿qué criterio utilizó para llegar a esa conclusión?

9. Realice los gráficos exploratorios de residuos para comprobación de supuestos del modelo

hist(TRanov$residuals)

boxplot(TasaR$TR~TasaR$DDC, main="TASA RESPIRATORIA", xlab="DDC", ylab="TR",col=4)# Box plot por tratamientos

boxplot(TRanov$residuals~TasaR$DDC,main="Residuales del modelo",col=4)

¿Cómo describiría el comportamiento de los residuos en el histograma? ¿Es una distribución normal? ¿Qué diferencia encuentra entre el boxplot por tratamientos y el de residuos del modelo? ¿Qué información puede inferir de estas graficas con respecto a los supuestos del modelo?

###Gráficos de validación de supuestos sobre residuales###
split.screen(c(2,2))

[1] 1 2 3 4

screen(1)

# gráfico de residuales vs valores ajustados
plot(TRanov$fitted.values,TRanov$residuals,main="Residuals vs. Fitted", pch=20)
abline(h=0,lty=2)

screen(2)
# gráfico de residuales vs niveles de tratamientos
plot(TasaR$DDC,TRanov$residuals,main="Residuals vs. Levels")   

screen(3)
# gráfico de residuales vs valores en el tiempo
plot(1:100,TRanov$residuals, main="Residuals vs. Time order",pch=20)  

abline(h=0,lty=2)
screen(4)

qqnorm(TRanov$residuals,pch=20)  
qqline(TRanov$residuals) # linea de probabilidad normal gráfico de residuales Q-Q

10. Compruebe los supuestos del modelo con las pruebas de normalidad y homocedasticidad.

### prueba shapiro prueba de normalidad ###
shapiro.test(TRanov$residuals)

    Shapiro-Wilk normality test

data:  TRanov$residuals
W = 0.97677, p-value = 0.0743

### prueba Kolmogorov-Smirnov prueba de normalidad ###
ks.test(TRanov$residuals, "pnorm", mean(TRanov$residuals), sd(TRanov$residuals))

    One-sample Kolmogorov-Smirnov test

data:  TRanov$residuals
D = 0.059213, p-value = 0.8745
alternative hypothesis: two-sided

### prueba bartlett prueba de homocedasticidad ###
bartlett.test(TRanov$residuals,TasaR$DDC)

    Bartlett test of homogeneity of variances

data:  TRanov$residuals and TasaR$DDC
Bartlett's K-squared = 7.8656, df = 9, p-value = 0.5477

¿Qué puede inferir de las p-valor de las pruebas? ¿Hay normalidad de los datos? ¿Hay homocedasticidad?

11. Realice la prueba de tukey para verificar cuales son les tratamientos (Días) presentan diferencias significativas con respecto a la variable de tasa respiratoria

HSD.test(TRanov,"TasaR$DDC",console = T)

Study: TRanov ~ "TasaR$DDC"

HSD Test for TasaR$TR 

Mean Square Error:  95.24518 

TasaR$DDC,  means

Alpha: 0.05 ; DF Error: 90 
Critical Value of Studentized Range: 4.588313 

Minimun Significant Difference: 14.16037 

Treatments with the same letter are not significantly different.

¿Hay diferencias significativas? ¿Entre cuales tratamientos? ¿Cómo este resultado se relaciona con el ANOVA realizado previamente?

12. Realice la gráfica Tasa Respiratoria vs. DDC que incluya las diferencias significativas.

TR=boxplot(TasaR$TR~TasaR$DDC,   main = "Tasa Respiratoria",
           ylab = "mg CO2 * kg-1 * h-1",
           xlab = "Días despues de cosecha", ylim=c(0,84),  col=brewer.pal(n = 8, name = "Set2"))
text(1,41,c("e"))

text(2,54,c("de"))
text(3,54,c("bcd"))

text(4,77,c("a"))
text(5,84,c("a"))

text(6,70,c("ab"))
text(7,70,c("bcd"))

text(8,65,c("bcd"))
text(9,63,c("bcd"))

text(10,63,c("cd"))

13. Ahora vamos a ver una variable diferente en la que los supuestos no se cumplen y es necesario realizar una transformación de los datos. Para eso organice los datos para la variable de Pérdida de Peso Fresco en columnas y llame el archivo de excel desde R.

Peso_Fresco<- read_excel("C:/Users/Juliana Miranda/Documents/Universidad/TESIS/RESULTADOS/R/Perdida_peso_fresco.xlsx")
Peso_Fresco

summary(Peso_Fresco)

      DDC          PPF              PPF1.2          PPF_BOX      
 Min.   : 4   Min.   : 0.7874   Min.   :0.3937   Min.   :0.9804  
 1st Qu.: 6   1st Qu.: 3.8819   1st Qu.:1.9409   1st Qu.:1.1188  
 Median : 8   Median : 6.7167   Median :3.3583   Median :1.1708  
 Mean   : 8   Mean   : 6.5402   Mean   :3.2701   Mean   :1.1494  
 3rd Qu.:10   3rd Qu.: 9.2884   3rd Qu.:4.6442   3rd Qu.:1.2027  
 Max.   :12   Max.   :13.6531   Max.   :6.8266   Max.   :1.2416  

14. Determine el peso fresco como un factor

### El tratamiento debe ser un factor ####
Peso_Fresco$DDC=as.factor(Peso_Fresco$DDC)

15. Cree el modelo y realice el análisis de varianza

PPFanov<-aov(Peso_Fresco$PPF~Peso_Fresco$DDC)
summary(PPFanov)

                Df Sum Sq Mean Sq F value Pr(>F)    
Peso_Fresco$DDC  8 1020.9  127.62   116.7 <2e-16 ***
Residuals       81   88.6    1.09                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

¿Qué se puede inferir de esta información? ¿qué criterio utilizó para llegar a esa conclusión?

16. Verifique los supuestos con las pruebas de normalidad y homocedasticidad

### prueba shapiro prueba de normalidad ###
shapiro.test(PPFanov$residuals)

    Shapiro-Wilk normality test

data:  PPFanov$residuals
W = 0.94696, p-value = 0.001091

### prueba Kolmogorov-Smirnov prueba de normalidad ###
ks.test(PPFanov$residuals, "pnorm", mean(PPFanov$residuals), sd(PPFanov$residuals))

Warning in ks.test(PPFanov$residuals, "pnorm", mean(PPFanov$residuals),  :
  ties should not be present for the Kolmogorov-Smirnov test

    One-sample Kolmogorov-Smirnov test

data:  PPFanov$residuals
D = 0.14315, p-value = 0.05001
alternative hypothesis: two-sided

### prueba de homogeneidad de varianzas ####
bartlett.test(PPFanov$residuals,Peso_Fresco$DDC)

    Bartlett test of homogeneity of variances

data:  PPFanov$residuals and Peso_Fresco$DDC
Bartlett's K-squared = 48.737, df = 8, p-value = 7.141e-08

¿Se cumplen los supuestos? ¿qué criterio utilizó para llegar a esa conclusión?

17. Intente transformar los datos con transformaciones comunes como la raíz cuadrada, logaritmo, etc. e intente comprobar nuevamente los supuestos. Para eso debe modificar el Excel creado y añadir una nueva columna con las variables transformadas, por ejemplo, a la 1/2 (la raíz cuadrada). Repita el proceso del punto 13 al 16 de este taller.

PPF1.2anov<-aov(Peso_Fresco$PPF1.2~Peso_Fresco$DDC)
summary(PPF1.2anov)

                Df Sum Sq Mean Sq F value Pr(>F)    
Peso_Fresco$DDC  8 255.23   31.90   116.7 <2e-16 ***
Residuals       81  22.14    0.27                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

### prueba shapiro prueba de normalidad ###
shapiro.test(PPF1.2anov$residuals)

    Shapiro-Wilk normality test

data:  PPF1.2anov$residuals
W = 0.94696, p-value = 0.001091

### prueba Kolmogorov-Smirnov prueba de normalidad ###
ks.test(PPF1.2anov$residuals, "pnorm", mean(PPF1.2anov$residuals), sd(PPF1.2anov$residuals))

Warning in ks.test(PPF1.2anov$residuals, "pnorm", mean(PPF1.2anov$residuals),  :
  ties should not be present for the Kolmogorov-Smirnov test

    One-sample Kolmogorov-Smirnov test

data:  PPF1.2anov$residuals
D = 0.14315, p-value = 0.05001
alternative hypothesis: two-sided

### prueba de homogeneidad de varianzas ####
bartlett.test(PPF1.2anov$residuals,Peso_Fresco$DDC)

    Bartlett test of homogeneity of variances

data:  PPF1.2anov$residuals and Peso_Fresco$DDC
Bartlett's K-squared = 48.737, df = 8, p-value = 7.141e-08

¿Mejoraron los valores de los supeustos? ¿Se cumplen?

18. Implemente el método Box-Cox para obtener el labda óptimo para la transformación de los datos.

datos<-data.frame(Peso_Fresco$PPF,Peso_Fresco$DDC)
bc1<-boxcox(Peso_Fresco$PPF~Peso_Fresco$DDC, lambda = seq(-2, 3, 1/5), data=datos, plotit = TRUE,
           eps = 1/50 , xlab = expression(lambda), ylab = "log-Likelihood")

summary(powerTransform(Peso_Fresco$PPF~Peso_Fresco$DDC))

bcPower Transformation to Normality 
   Est Power Rounded Pwr Wald Lwr Bnd Wald Upr Bnd
Y1    0.0828           0      -0.1059       0.2716

Likelihood ratio test that transformation parameter is equal to 0
 (log transformation)

Likelihood ratio test that no transformation is needed

¿Qué muestra la gráfica? ¿Qué valor lambda se obtiene?

19. Realice la transformación de los datos con su nuevo valor lambda. Nuevamente modifique el excel original y vuelva a cargar los datos

##Datos transformados ^(0.0828)######

PPFboxanov<-aov(Peso_Fresco$PPF_BOX~Peso_Fresco$DDC)
summary(PPFboxanov)

                Df Sum Sq Mean Sq F value Pr(>F)    
Peso_Fresco$DDC  8 0.4205 0.05257   268.3 <2e-16 ***
Residuals       81 0.0159 0.00020                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

### prueba shapiro prueba de normalidad ###
shapiro.test(PPFboxanov$residuals)

    Shapiro-Wilk normality test

data:  PPFboxanov$residuals
W = 0.95825, p-value = 0.0056

### prueba Kolmogorov-Smirnov prueba de normalidad ###
ks.test(PPFboxanov$residuals, "pnorm", mean(PPFboxanov$residuals), sd(PPFboxanov$residuals))

Warning in ks.test(PPFboxanov$residuals, "pnorm", mean(PPFboxanov$residuals),  :
  ties should not be present for the Kolmogorov-Smirnov test

    One-sample Kolmogorov-Smirnov test

data:  PPFboxanov$residuals
D = 0.088869, p-value = 0.4759
alternative hypothesis: two-sided

### prueba de homogeneidad de varianzas ####
bartlett.test(PPFboxanov$residuals,Peso_Fresco$DDC)

    Bartlett test of homogeneity of variances

data:  PPFboxanov$residuals and Peso_Fresco$DDC
Bartlett's K-squared = 1.0237, df = 8, p-value = 0.9981

¿El valor lambda permitió corregir los supuestos del modelo? ¿Con qué criterio llegó a esa conclusión?

20. Con los supuestos cumplidos, realice el análisis de varianza con los datos transformados

HSD.test(PPFboxanov,"Peso_Fresco$DDC",console = T)

Study: PPFboxanov ~ "Peso_Fresco$DDC"

HSD Test for Peso_Fresco$PPF_BOX 

Mean Square Error:  0.0001959574 

Peso_Fresco$DDC,  means

Alpha: 0.05 ; DF Error: 81 
Critical Value of Studentized Range: 4.507391 

Minimun Significant Difference: 0.0199529 

Treatments with the same letter are not significantly different.

21. Grafique los resultados

PPFg=boxplot(Peso_Fresco$PPF~Peso_Fresco$DDC,   main = "Pérdida peso fresco",
          ylab = "% Pérdida peso fresco",
          xlab = "Días despues de cosecha", col=brewer.pal(n = 8, name = "Set2"))
text(1,2,c("g"))

text(2,4,c("f"))
text(3,6,c("e"))

text(4,7.4,c("d"))
text(5,9,c("c"))

text(6,11,c("bc"))
text(7,12,c("ab"))

text(8,13,c("a"))
text(9,14,c("a"))

DCA con Submuestreo

22. Para este análisis se usará la variable de Firmeza. Nuevamente organice las variables de forma que queden en columnas de DDC, Firmeza y muestra, esta última en este caso corresponde a los valores del submuestreo (medidas hechas al mismo fruto). Una vez creado el Excel, llame el documento desde R.

Textura<- read_excel("C:/Users/Juliana Miranda/Documents/Universidad/TESIS/RESULTADOS/R/Txt.xlsx")
Textura

summary(Textura)

     DDC                 FIRM           MUESTRA 
 Length:36          Min.   : 4.426   Min.   :1  
 Class :character   1st Qu.: 6.372   1st Qu.:1  
 Mode  :character   Median : 8.220   Median :2  
                    Mean   :16.912   Mean   :2  
                    3rd Qu.:33.708   3rd Qu.:3  
                    Max.   :41.948   Max.   :3  

23.Establezca tanto los estados poscosecha como las medidas del submuestreo como factores

ESTADO<- as.factor(Textura$DDC)
MUESTRA<- as.factor(Textura$MUESTRA) 

24. Cree el modelo estadístico para un DCA con submuestreo

FIRMdca2<-aov(Textura$FIRM~ESTADO+MUESTRA%in%ESTADO)
summary(FIRMdca2)

               Df Sum Sq Mean Sq F value Pr(>F)    
ESTADO          2   6882    3441 843.590 <2e-16 ***
ESTADO:MUESTRA  6     30       5   1.211  0.331    
Residuals      27    110       4                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

25. Cree los residuos del modelo y determine si se cumplen los supuestos

resFIRM<-FIRMdca2$residuals #RESIDUOS 
resFIRM

         1          2          3          4          5          6          7          8          9         10         11 
-1.3815000 -1.3135000  0.7305000  1.9645000 -4.5212500 -3.1202500  6.0197500  1.6217500  0.6710000  3.1380000 -2.2930000 
        12         13         14         15         16         17         18         19         20         21         22 
-1.5160000  1.0058250 -0.8537750  0.6624250 -0.8144750 -0.6815500  0.2950500 -0.7906500  1.1771500  1.2299775 -0.4369225 
        23         24         25         26         27         28         29         30         31         32         33 
 0.7340775 -1.5271325  0.7579500 -0.0030500 -1.3038500  0.5489500 -0.2633250 -0.2633250  0.3147750  0.2118750  0.2222750 
        34         35         36 
 0.9556750  0.0049750 -1.1829250 

shapiro.test(resFIRM)

    Shapiro-Wilk normality test

data:  resFIRM
W = 0.92979, p-value = 0.02467

bartlett.test(resFIRM,ESTADO)

    Bartlett test of homogeneity of variances

data:  resFIRM and ESTADO
Bartlett's K-squared = 24.347, df = 2, p-value = 5.166e-06

ks.test(FIRMdca2$residuals, "pnorm", mean(FIRMdca2$residuals), sd(FIRMdca2$residuals))

Warning in ks.test(FIRMdca2$residuals, "pnorm", mean(FIRMdca2$residuals),  :
  ties should not be present for the Kolmogorov-Smirnov test

    One-sample Kolmogorov-Smirnov test

data:  FIRMdca2$residuals
D = 0.13292, p-value = 0.5482
alternative hypothesis: two-sided

26. Continue el análisis determinando la transformación apropiada para el cual se cumplan los supuestos. Determine las diferencias signficativas entre los estados poscosecha. Grafique los resultados.