rm(list=ls())
library(ggplot2)
library(rgl)
library(lattice)
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, union
library(MASS)
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
library(readxl)
library(pander)
library(shape)
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
## 
## Attaching package: 'xts'
## The following objects are masked from 'package:dplyr':
## 
##     first, last
## 
## Attaching package: 'PerformanceAnalytics'
## The following object is masked from 'package:graphics':
## 
##     legend
library(tidyverse)
## -- Attaching packages ----------------------------------------- tidyverse 1.3.0 --
## v tibble  2.1.3     v purrr   0.3.3
## v tidyr   1.0.2     v stringr 1.4.0
## v readr   1.3.1     v forcats 0.5.0
## -- Conflicts -------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x xts::first()    masks dplyr::first()
## x dplyr::lag()    masks stats::lag()
## x xts::last()     masks dplyr::last()
## x MASS::select()  masks dplyr::select()
library(survival)
library(coin)
library(pastecs)
## 
## Attaching package: 'pastecs'
## The following object is masked from 'package:tidyr':
## 
##     extract
## The following objects are masked from 'package:xts':
## 
##     first, last
## The following objects are masked from 'package:dplyr':
## 
##     first, last
library(symbols)
library(scatterplot3d)
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:purrr':
## 
##     some
## The following object is masked from 'package:dplyr':
## 
##     recode
library(outliers)

Ejercicio 1

Se midió la conductancia estomática (gs: mol/m2s) en dos cultivares de papa diploide (Colombia y Ocarina) bajo una condición de déficit de riego.

Determinar al 95% de nivel de confianza si las dos medias obtenidas para los cultivares son estadísticamente iguales. Utilice la información del artículo mostrado en clase para decidir si las varianzas pueden considerarse iguales o no.

colom = (c(0.45, 0.41,0.47, 0.46,0.39,0.44,0.48, 0.42, 0.44, 0.48, 0.50, 0.47,0.44,0.52))

ocari= (c(0.28,0.25,0.32,0.34,0.36,0.40,0.37,0.36,0.39,0.41,0.37,0.42,0.41))

*Se debe probar los supuestosde normalidad en los datos

par(mfrow=c(1,2))

histogram(colom)

histogram(ocari, na.rm = T)

par(mfrow=c(1,2))
qqnorm(colom)
qqline(colom, col="red")
qqnorm(ocari)
qqline(ocari, col="red")

Los gráficos no son suficientes para concluir que los datos tiene distribución normal.

shapiro.test (colom)
## 
##  Shapiro-Wilk normality test
## 
## data:  colom
## W = 0.98666, p-value = 0.9971

Cómo el resutado del p-valor es de 0.91 , una cifra mucho mayor a 0,05, podemos decir que no hay evidencia para rechazar la hipotesis nula que en este caso es que los datos tienen distribucion normal.

shapiro.test (ocari)
## 
##  Shapiro-Wilk normality test
## 
## data:  ocari
## W = 0.90993, p-value = 0.183

Cómo el resutado del p-valor es de 0.183 , una cifra mayor a 0,05, podemos decir que no hay evidencia para rechazar la hipotesis nula que en este caso es que los datos tienen distribucion normal.

var_colom=var(colom)
var_ocari=var((ocari))
var_colom
## [1] 0.001242308
var_ocari
## [1] 0.00265
vp1= var.test(x= colom, y= ocari , alternative= "t",conf.level = 0.95)
vp1
## 
##  F test to compare two variances
## 
## data:  colom and ocari
## F = 0.4688, num df = 13, denom df = 12, p-value = 0.1899
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.1447228 1.4781939
## sample estimates:
## ratio of variances 
##          0.4687954
ifelse(vp1$p.value<0.05, 
       'Rechazo Ho: Las varianzas no son iguales',
        'No rechazo Ho: Las varianzas son iguales')
## [1] "No rechazo Ho: Las varianzas son iguales"

Según el artículo las varianzas no se pueden considerar iguales porque una es mas del doble de la otra.

Hipótesis

\[H_0: \mu_{po}=\mu_{pc}\] \[H_a: \mu_{po}≠ \mu_{pc}\]

prb_1 = t.test(colom,ocari, alternative = "t", mu=0, conf.level = 0.95,var.equal = F)
prb_1
## 
##  Welch Two Sample t-test
## 
## data:  colom and ocari
## t = 5.5539, df = 21.041, p-value = 1.63e-05
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.05943233 0.13056767
## sample estimates:
## mean of x mean of y 
##     0.455     0.360
prb_1$p.value
## [1] 1.629826e-05
ifelse(prb_1$p.value<0.05, 
       'Rechazo Ho: Las medias no son iguales',
         'No rechazo Ho: Las medias no son iguales')
## [1] "Rechazo Ho: Las medias no son iguales"

**Conclusion Al correr la prueba T student el pvalor es por mucho menor de 0,05, por lo que tenemos evidencia estadistica para rechazar la hipotesis nula,es decir las medias de la conductancia estomatica en las variedades de papa diploide expuestas a un misma condicion de riego insuficiente no son iguales. Esto nos ayuda a concluir que variedad Colombia tiene menos limitaciones ante una posible temporadas de sequia y tiene una mejor produccion de fotoasimilados y por ende un mejor rendimiento que la la variedad Ocarina.

Ejercicio 2

Se propuso un plan de fertilización en papa criolla y se midió a los 45 y 77 días después de la siembra el peso de tubérculos (Kg/ha) .

dia_45=(c(60,64,65,66,66,66,67,67,67,68,68,68,68,69,69,72))
dia_77=(c(759,832,843,850,840,920,875,910,800,834,790,840,
          812,873,905,832))

Muestras dependientes

En este experimento tenemos unos datos pareados. Vamos a realizar un histograma de las diferencias entre la primer y segunda medición.

daf <- data.frame (dia_45, dia_77)

Debemos determinar al 95% de nivel de confianza si se incrementó la medida de rendimiento en las dos evaluaciones registradas.

Debemos revisar los supuestos de las pruebas.

shapiro.test(dia_77)
## 
##  Shapiro-Wilk normality test
## 
## data:  dia_77
## W = 0.9645, p-value = 0.7435
shapiro.test(dia_45)
## 
##  Shapiro-Wilk normality test
## 
## data:  dia_45
## W = 0.91678, p-value = 0.1497

Los datos presentan distribucion normal, y procedemos a aplicar la prueba T Student

Hipótesis

\[H_0: \mu_{d77}>\mu_{d45}\]

\[H_a: \mu_{d77}≤ \mu_{d45}\]

prueba_2 = t.test(x = dia_45, y = dia_77, alternative = 't', mu = 0, conf.level = 0.95)
prueba_2
## 
##  Welch Two Sample t-test
## 
## data:  dia_45 and dia_77
## t = -70.387, df = 15.105, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -801.3519 -754.2731
## sample estimates:
## mean of x mean of y 
##   66.8750  844.6875
ifelse(prueba_2$p.value<0.05,
      'Rechazo Ho: El rendimiento del dia77 es diferente al rendimiento del dia45',
      'No rechazo Ho: El rendimiento del dia77 es igual al rendimiento del dia45')
## [1] "Rechazo Ho: El rendimiento del dia77 es diferente al rendimiento del dia45"

El p.valor es menor que 0,05 ,por lo que se rechaza la hipotesis nula.

cor1= cor.test(dia_45, dia_77)
cor1
## 
##  Pearson's product-moment correlation
## 
## data:  dia_45 and dia_77
## t = 1.3274, df = 14, p-value = 0.2056
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1934072  0.7120442
## sample estimates:
##       cor 
## 0.3343536

El valor de la hallado al correr la test de pearson es de 0,33. Es un valor pequeño que nos indica que no existe una relacion directa, positiva debil entre las dos mediciones y este resutado tiene sentido ya que en el la acumulación de azucares en los tuberculos hay muchos factores que pueden afectar.

Coeficiente de determinación

0.3343^2
## [1] 0.1117565

Solo se explican el 11% de los datos con esta relación.

chart.Correlation(daf)

Con este gráfico podemos observar el comportamiento de las mediciones, no existe una relación lineal.

###Cambio relativo porcentual promedio

colMeans(daf)
##   dia_45   dia_77 
##  66.8750 844.6875

Calculamos las medias de cada medición para hallar el cambio relativo.

((mean(daf$dia_77)-mean(daf$dia_45))/(mean(daf$dia_77))*100)
## [1] 92.08287

El cambio relativo promedio de las medias de las mediciones es de 92.08%, lo que significa que en el término de 32 dias, el incremento de peso por hectarea ha incrementado elevado.

**Conclusion En este ejercicio se puede concluir que la media del peso seco en el dia 77 no es el mismo que el dia 45 (se rechaza la hipotesis nula), por medio del coeficiente de Pearson sabemos estas muestras estan debilmente relacionadas, es decir elcambio en el rendimiento no esta ligado al primer rendimiento sino hay mucho otros factores responsables de esa variación (de los que no tenemos información).

Ejercicio 3 .Prueba de Wilcoxon de la suma de rangos-Dos muestras independientes

Se está evaluando la calidad de frito mediante la textura de las hojuelas de papa criolla en dos tipos de aceite (palma y maíz) utilizado para freír en condiciones controladas de tiempo y temperatura. Al final se recolectaron las hojuelas y se evaluó en una escala diagramática la calidad de frito (escala de 1 a 5, desde (1) no crujiente hasta (5) bastante crujientes).

En este caso tenemos una muestra pequeña, para determinar si tiene distrubución normal. Se usará una prueba no paramétrica, al tratarse de dato cualitativos dispuestos en rangos.

palma=c(3,4,3,4,4,3,3,4,4,3,4,4,2,4,3,4,3,3,3,4,4)
maiz= c(3,4,4,4,4,4,3,4,3,4,4,4,4,3,4,4,4,3,3,4,3)
summary(palma)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   2.000   3.000   4.000   3.476   4.000   4.000
summary(maiz)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   3.000   3.000   4.000   3.667   4.000   4.000

\[H_0: \median_{oil}=0\]

\[H_a: \median_{oil} ≠0\] Se aplican pruebas para comprobar suspuestos .

fligner.test(x = list(palma,maiz))
## 
##  Fligner-Killeen test of homogeneity of variances
## 
## data:  list(palma, maiz)
## Fligner-Killeen:med chi-squared = 1.3429, df = 1, p-value = 0.2465

No hay evidencia para negar la normalidad de las varianzas

n1 <- length(palma)
n2 <- length(maiz)
par(mfrow=c(1,2))
plot(palma)
plot(maiz)

aceite= gl(2,21,42, c("Palma", "Maiz"))
datos=data.frame(c(palma,maiz),aceite)
colnames(datos)=c("Cal_fritura", "aceite")
ggplot(datos, aes(x = Cal_fritura , fill = aceite)) +
  geom_density(alpha = 0.6)

En l gráfica se puede observar las medidas de la calidad de la fritura con los dos tipos de aceite.

wilcox.test(x = palma, y = maiz, alternative = "two.sided", mu = 0,
            paired = F, conf.int = 0.95, correct = F, exact = F)
## 
##  Wilcoxon rank sum test
## 
## data:  palma and maiz
## W = 185.5, p-value = 0.3042
## alternative hypothesis: true location shift is not equal to 0
## 95 percent confidence interval:
##  -9.999919e-01  1.939953e-05
## sample estimates:
## difference in location 
##          -8.889014e-05
datos <- data.frame(grupo = rep(c("A", "B"), c(21, 21)),
                    valores = c(palma,maiz))
wilcox_test(valores ~ grupo, data = datos, distribution = "exact", conf.int=0.95)
## 
##  Exact Wilcoxon-Mann-Whitney Test
## 
## data:  valores by grupo (A, B)
## Z = -1.0275, p-value = 0.3976
## alternative hypothesis: true mu is not equal to 0
## 95 percent confidence interval:
##  0 0
## sample estimates:
## difference in location 
##                      0
tamanyo_efecto <- 1.0275/sqrt(n1 + n2)
tamanyo_efecto
## [1] 0.1585467

**Conclusion En la litertura encontramos dos maneras de correr la prueba de Wilconson, las dos nos arrojan un p.valor mayor a 0.05,por lo que no hay evidencia estadistica para rechazar hipotesis nula, es decir las diferencia de las medianas de las muestras es igual a cero, por lo que se puede decir que la calidad de la fritura de las hojuelas de papa no dependen del aceite en el que se fria, sino que estas diferencias pueden radicar en las carasteristicas del tuberculo.

Ejercicio 4

Prueba de Wilcoxon de la suma de rangos-Dos muestras pareadas Suponga que del ejercicio anterior se seleccionó el aceite de maíz y se desarrolló un segundo experimento para controlar la temperatura de almacenamiento de papa criollo. Las temperaturas de almacenamiento fueron 4 y 12 °C y se utilizó la escala de color CIELab*.

eje_4 <- read_excel("c:/Users/Usuario/Desktop/acei4_maiz.xlsx")
eje_12 <- read_excel("c:/Users/Usuario/Desktop/ace12_maiz.xlsx")
L= c(69.26,68.15,69.17,68.88,70.01,70.15,70.66,68.68,71,72.18,69.15,70,68.64,68.12,68.12,62.2,60.45,63.12,61.64,61.25,62.55,64.12,65.65,66.87,65.11,66.14,62.64,61.97,60.58,60.68)
a= c(-1.31,-1.25,-1.47,-1.35,-1.32,-1.15,-1.25,-1.29,-1.42,-1.45,-1.29,-1.22,-1.19,-1.25,-1.25,0.81,0.78,0.55,0.81,0.77,0.69,0.59,0.55,0.42,0.39,0.41,0.37,0.35,0.34,0.34)
b= c(28.68,27.66,28.02,27.66,27.66,26.88,26.25,26.26,28.15,30,28.24,25.59,24.69,25.56,26.26,37.31,35.9,36.36,36.12,36.45,35.99,36.14,36.14,35.55,34.77,32.32,31.96,30.17,36.65,37.15)
temp= gl(2,15,30,labels = c("4°C", "12°C"))
datos3=data.frame(L,a,b,temp)
par(mfrow= c(3, 3), cex = 0.6, mar = c(3,3,3,1))
ggplot(datos3, aes(x = L, fill = temp)) +
  geom_density(alpha = 0.6)

ggplot(datos3, aes(x = a, fill = temp)) +
  geom_density(alpha = 0.6)

ggplot(datos3, aes(x = b, fill = temp)) +
  geom_density(alpha = 0.6)

Usaremos eltest de Wilcolson

wilcox.test(x =eje_4$L, y = eje_12$L, alternative = "two.sided", mu = 0,
            paired = T, conf.int = 0.98, correct = F, exact = F)
## 
##  Wilcoxon signed rank test
## 
## data:  eje_4$L and eje_12$L
## V = 120, p-value = 0.000655
## alternative hypothesis: true location shift is not equal to 0
## 95 percent confidence interval:
##  5.304912 7.420047
## sample estimates:
## (pseudo)median 
##       6.883541
wilcox.test(x =eje_4$a, y = eje_12$a, alternative = "two.sided", mu = 0,
            paired = T, conf.int = 0.98, correct = F, exact = F)
## 
##  Wilcoxon signed rank test
## 
## data:  eje_4$a and eje_12$a
## V = 0, p-value = 0.000637
## alternative hypothesis: true location shift is not equal to 0
## 95 percent confidence interval:
##  -1.965046 -1.714983
## sample estimates:
## (pseudo)median 
##      -1.840011
wilcox.test(x =eje_4$b, y = eje_12$b, alternative = "two.sided", mu = 0,
            paired = T, conf.int = 0.98, correct = F, exact = F)
## 
##  Wilcoxon signed rank test
## 
## data:  eje_4$b and eje_12$b
## V = 0, p-value = 0.000655
## alternative hypothesis: true location shift is not equal to 0
## 95 percent confidence interval:
##  -9.335069 -6.884976
## sample estimates:
## (pseudo)median 
##      -8.287353

Vemos que al aplicar la prueba de wilcoxon a los componentes de las escala de color en las temperaturas 4 y 12 grados se halla que el L* y b* tienen valores iguales (0,000655), miientras que en a* el valor es de 0,000637. Por lo que se puede decir que el color de los chips de papa variara levemente segun la temperatura de fritura.

E =sqrt((eje_12$L-eje_4$L)^2)+(eje_12$a-eje_4$a)^2+(eje_12$b-eje_4$b)
E
##  [1] 20.1844 20.0609 18.4704 20.3656 21.9181 20.0956 19.8156 16.2956 14.9156
## [10] 15.2256  9.9800 16.2581 14.5216 21.1581 20.8581
wilcox.test(E, mu = 0, alternative = "two.sided",conf.level = 0.95,
            correct= F, exact= F)
## 
##  Wilcoxon signed rank test
## 
## data:  E
## V = 120, p-value = 0.000655
## alternative hypothesis: true location is not equal to 0

El p.valor resultante de la prueba Wicoxon realizada a los valores transformados es igual a la de la prueba en “L” y en “b”, pero difiere en “a”. En la escalas de colores “a” corresponde a las tonalidades entre el rojo y el azul.

x4= eje_4$L
y4= eje_4$a
z4= eje_4$b
x12= eje_12$L
y12= eje_12$a
z12= eje_12$b
par(mfrow=c(1,2))
scatterplot3d(x4,y4,z4,pch= 15, "pink")
scatterplot3d(x12,y12,z12,pch= 16, "lightblue")

**Conclusion

Al analizar los datos de las escalas de colores con las dos temperturas experimentadas hallamos que hay una variación pequeña en la escala “a *”, mientras que en “L*” y “b*” elresultado es el mismo, por lo tanto concluir que al variar la temperatura de fritura, habra un muy leve cambio en el color del chip,ya que al realizar las transformaciones obtenemos el mismo valor que en “L*” y "b*".

Ejercicio 5

Prueba F- AOV-FSCA-B

Un laboratorio emplea 3 métodos (Bray,Olsen; Mehlich-3) para determinar el contenido de fósforo en suelos. Surge la pregunta: “¿Difieren las medias determinaciones entre los métodos de análisis?”

fosforo <- read_excel("c:/Users/Usuario/Desktop/soil.xlsx")
Bray =c(7.1,6.8,6.6,6.7,6.8,6.7,6.9,6.8,6.7,6.6,6.9,
       7.4,6.6,6.7,6.8,7,6.7,6.7,7.1,6.6,6.6,6.9)

Olsen =c(6,6.2,6.4,6.6,6.9,6.6,6.4,6.4,6.5,6.3,6.4,
        6.6,6.6,6.7,6.4,6.2,6.6,6.3,5.8,6.5,6.2,6.4)

Mehlich =c(9.1,7.1,7.8,7.3,7.6,7.8,7.4,7.3,7.3,7.1,8,
          8.7,7.6,7.8,7.8,7.2,7.2,7.7,7.3,7.5,7.1,7.1)
metod= gl(3,22,66, c("Bray", "Olsen", "Mehlich"))
dato5=data.frame(c(Bray, Olsen, Mehlich), metod)
dat5=data.frame(f=fosforo, m=metod)
colnames(dato5)=c("Fosforo","Metodo")
ggplot(dato5, aes(x = Fosforo, fill = metod)) +
  geom_density(alpha = 0.6)

stat.desc(fosforo)
##                      Bray        Olsen      Mehlich
## nbr.val       22.00000000  22.00000000  22.00000000
## nbr.null       0.00000000   0.00000000   0.00000000
## nbr.na         0.00000000   0.00000000   0.00000000
## min            6.60000000   5.80000000   7.10000000
## max            7.40000000   6.90000000   9.10000000
## range          0.80000000   1.10000000   2.00000000
## sum          149.70000000 141.00000000 166.80000000
## median         6.75000000   6.40000000   7.45000000
## mean           6.80454545   6.40909091   7.58181818
## SE.mean        0.04338356   0.05134936   0.10923511
## CI.mean.0.95   0.09022105   0.10678685   0.22716686
## var            0.04140693   0.05800866   0.26251082
## std.dev        0.20348692   0.24084987   0.51235810
## coef.var       0.02990456   0.03757941   0.06757721
$$H_0: \tau_{Bray}\=\tau_{Olsen}\=\tau_{Mehlib}\=0$$

\[H_a:H_0 es falsa\]

anova <- aov(dato5$Fosforo ~ dato5$Metodo)
summary(anova)
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## dato5$Metodo  2  15.66   7.831   64.91 4.97e-16 ***
## Residuals    63   7.60   0.121                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

F value es un valor muy elevado, lo que nos indica que se rechaza la hipotesis nula.

res1 = anova$residuals
res1
##            1            2            3            4            5            6 
##  0.295454545 -0.004545455 -0.204545455 -0.104545455 -0.004545455 -0.104545455 
##            7            8            9           10           11           12 
##  0.095454545 -0.004545455 -0.104545455 -0.204545455  0.095454545  0.595454545 
##           13           14           15           16           17           18 
## -0.204545455 -0.104545455 -0.004545455  0.195454545 -0.104545455 -0.104545455 
##           19           20           21           22           23           24 
##  0.295454545 -0.204545455 -0.204545455  0.095454545 -0.409090909 -0.209090909 
##           25           26           27           28           29           30 
## -0.009090909  0.190909091  0.490909091  0.190909091 -0.009090909 -0.009090909 
##           31           32           33           34           35           36 
##  0.090909091 -0.109090909 -0.009090909  0.190909091  0.190909091  0.290909091 
##           37           38           39           40           41           42 
## -0.009090909 -0.209090909  0.190909091 -0.109090909 -0.609090909  0.090909091 
##           43           44           45           46           47           48 
## -0.209090909 -0.009090909  1.518181818 -0.481818182  0.218181818 -0.281818182 
##           49           50           51           52           53           54 
##  0.018181818  0.218181818 -0.181818182 -0.281818182 -0.281818182 -0.481818182 
##           55           56           57           58           59           60 
##  0.418181818  1.118181818  0.018181818  0.218181818  0.218181818 -0.381818182 
##           61           62           63           64           65           66 
## -0.381818182  0.118181818 -0.281818182 -0.081818182 -0.481818182 -0.481818182

Aplicamos las pruebas para validar supuestos

norm = shapiro.test(res1)
norm
## 
##  Shapiro-Wilk normality test
## 
## data:  res1
## W = 0.86946, p-value = 5.098e-06
home = bartlett.test(res1, dato5$Metodo)
home
## 
##  Bartlett test of homogeneity of variances
## 
## data:  res1 and dato5$Metodo
## Bartlett's K-squared = 21.061, df = 2, p-value = 2.67e-05

El p.valor en las dos pruebas es menor de 0.05 por lo que no se cumplen los supuestos.

*Se realizaran transformaciones para buscar la normalidad en los datos.

**Transformación logaritmica

l_fos = log(x = dato5$Fosforo, base = 10)
anova1l = aov(l_fos ~ dato5$Metodo)

norml = shapiro.test(anova$residuals)
norml
## 
##  Shapiro-Wilk normality test
## 
## data:  anova$residuals
## W = 0.86946, p-value = 5.098e-06

Con la transformación logaritmica aun el p. valor es menor de 0,05, no se alcanza normalidad.

**Tranformacion raiz cuadrada

r_fos = sqrt(dato5$Fosforo)
anova1r = aov(r_fos ~dato5$Metodo)
normr = shapiro.test(anova1r$residuals)
normr
## 
##  Shapiro-Wilk normality test
## 
## data:  anova1r$residuals
## W = 0.89251, p-value = 3.276e-05

El p.valor es menor a 0,05 .Aun despues de la transformación , no hay normalidad.

**Transformación inversa

i_fosf = 1/dato5$Fosforo
anova1i = aov(i_fosf~dato5$Metodo)
normi = shapiro.test(anova1i$residuals)
normi
## 
##  Shapiro-Wilk normality test
## 
## data:  anova1i$residuals
## W = 0.9462, p-value = 0.006353

La tranformación logra normalidad.

Transformación Box-cox 
lb = powerTransform(dato5$Fosforo)
lb 
## Estimated transformation parameter 
## dato5$Fosforo 
##     -3.090281

Transformación de los datos

fos_t1 = (dato5$Fosforo)**lb$lambda
mod1_t1 = aov(fos_t1~dato5$Metodo)
normt = shapiro.test(mod1_t1$residuals)
normt
## 
##  Shapiro-Wilk normality test
## 
## data:  mod1_t1$residuals
## W = 0.96164, p-value = 0.03932

No se cumple normalidad.

En este ejemplo la transformción no nos ayuda a establecer la normalidad en los datos, por lo que se sospecha de estar ante datos atipicos. Se aplica el test de grubbs.

###Grubbs

d5 = grubbs.test(x = anova$residuals,two.sided = TRUE)
d5
## 
##  Grubbs test for one outlier
## 
## data:  anova$residuals
## G.45 = 4.43977, U = 0.69208, p-value = 8.785e-05
## alternative hypothesis: highest value 1.51818181818182 is an outlier

P valor < 0.05, rechazo \(H_0\), el valor 1,518 es un dato atipico.

Buscando el lugar del dato atipico en los datos.

which.max(anova$residuals)
## 45 
## 45
medias = tapply(dato5$Fosforo, dato5$Metodo, mean)
medias
##     Bray    Olsen  Mehlich 
## 6.804545 6.409091 7.581818
dato5$Fosforo[45] = medias[8]
anovac = aov(dato5$Fosforo~dato5$Metodo)
norc = shapiro.test(anovac$residuals)
summary (anovac)
##              Df Sum Sq Mean Sq F value Pr(>F)    
## dato5$Metodo  2 13.304   6.652   79.53 <2e-16 ***
## Residuals    62  5.186   0.084                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 1 observation deleted due to missingness
norc
## 
##  Shapiro-Wilk normality test
## 
## data:  anovac$residuals
## W = 0.92717, p-value = 0.0009079

Aun cambiando este dato atipico no se alcanza la normalidad. vamos a realizar las transformaciones.

**Segunda transformación Logaritmica

l2_fos = log(x = dato5$Fosforo, base = 10)
anoval2 = aov(l2_fos ~ dato5$Metodo)

norml2 = shapiro.test(anoval2$residuals)
norml2
## 
##  Shapiro-Wilk normality test
## 
## data:  anoval2$residuals
## W = 0.95626, p-value = 0.02197

No se alcanza normalidad

**Tranformacion raiz cuadrada

r2_fos = sqrt(dato5$Fosforo)
anovar2 = aov(r2_fos ~dato5$Metodo)
norr2 = shapiro.test(anovar2$residuals)
norr2
## 
##  Shapiro-Wilk normality test
## 
## data:  anovar2$residuals
## W = 0.94317, p-value = 0.004918

no se alcanza la normalidad

**Transformación inversa

i2_fosf = 1/dato5$Fosforo
anovai2 = aov(i2_fosf~dato5$Metodo)
nori2 = shapiro.test(anovai2$residuals)
nori2
## 
##  Shapiro-Wilk normality test
## 
## data:  anovai2$residuals
## W = 0.96948, p-value = 0.1083

Con la transformación inversa despues de haber cambiado el dato atipico por la media de ese metodo se alcanza la normalidad de los datos.

##Prueba de Tukey

tkf = TukeyHSD(anovac)
tkf
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = dato5$Fosforo ~ dato5$Metodo)
## 
## $`dato5$Metodo`
##                     diff        lwr        upr    p adj
## Olsen-Bray    -0.3954545 -0.6048443 -0.1860648 7.92e-05
## Mehlich-Bray   0.7049784  0.4931106  0.9168462 0.00e+00
## Mehlich-Olsen  1.1004329  0.8885651  1.3123007 0.00e+00

Esta prueba compara todos los pares posibles de los tratamientos, si tenemos un p <0.05 indica que los tratamientos comparados son diferentes. En el ejercicio vemos que no hay valores cercanos a 0.05 por lo que podemos decir que todas las pruebas son diferentes

qqnorm(dato5$Fosforo)
qqline(dato5$Fosforo, col= "violet")

xyplot(dato5$Fosforo~dato5$Metodo,main= "Residuales vs Nivelesde factores", xlab = "Combinaciones de nivel de factor", ylab = "Residuales Estandarizados")

boxplot(dato5$Fosforo~dato5$Metodo, xlab = NULL, ylab = "Contenido de fosforo")
points(medias, pch = 19, col = 'red')
abline(h = mean(dato5$Fosforo), lty = 2)

En el boxplot observamos la existencia de puntos outlayers,en todos los metodos.

data.aov <- aov (Fosforo ~ Metodo, dato5)
plot(data.aov)

Este código permite visualizar graficos de los residuales de los datos

anova(data.aov)
## Analysis of Variance Table
## 
## Response: Fosforo
##           Df  Sum Sq Mean Sq F value    Pr(>F)    
## Metodo     2 13.3040  6.6520  79.529 < 2.2e-16 ***
## Residuals 62  5.1858  0.0836                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1