Ejercicios diseño
Jcastillab
5/4/2021
CE_col=c(0.45,0.41,0.49,0.46,0.39,0.44,0.48,0.42,0.44,0.48,0.50,0.47,0.44,0.52)
CE_oca=c(0.28,0.25,0.32,0.34,0.36,0.40,0.39,0.36,0.39,0.41,0.37,0.42,0.41,NA)
dt=data.frame(CE_col,CE_oca);dt## CE_col CE_oca
## 1 0.45 0.28
## 2 0.41 0.25
## 3 0.49 0.32
## 4 0.46 0.34
## 5 0.39 0.36
## 6 0.44 0.40
## 7 0.48 0.39
## 8 0.42 0.36
## 9 0.44 0.39
## 10 0.48 0.41
## 11 0.50 0.37
## 12 0.47 0.42
## 13 0.44 0.41
## 14 0.52 NA\[H_0: \sigma^2_{Colombia} = \sigma^2_{Ocarina}\\ Ha:\sigma^2_{Colombia}\neq\sigma^2_{Ocarina}\]
prueba1=var.test(dt$CE_col,dt$CE_oca,ratio=1,alternative="t", conf.level=0.95)
ifelse(prueba1$p.value<0.05,"Rechazo Ho","No rechazo Ho")## [1] "No rechazo Ho"prueba2=t.test(dt$CE_col,dt$CE_oca,alternative = "t",mu = 0,var.equal = T,conf.level = 0.95)
ifelse(prueba2$p.value<0.05,"Rechazo Ho","No rechazo Ho")## [1] "Rechazo Ho"Varianzas iguales, medias diferentes
Ejercicio 2
Peso_45d=c(69,66,72,68,65,66,67,68,69,69,66,68,64,67,60,68)
Peso_77d=c(873,850,832,834,843,840,895,790,905,910,920,840,832,800,759,812)
rend=data.frame(Peso_45d,Peso_77d); rend## Peso_45d Peso_77d
## 1 69 873
## 2 66 850
## 3 72 832
## 4 68 834
## 5 65 843
## 6 66 840
## 7 67 895
## 8 68 790
## 9 69 905
## 10 69 910
## 11 66 920
## 12 68 840
## 13 64 832
## 14 67 800
## 15 60 759
## 16 68 812media_45d = mean(Peso_45d); media_45d## [1] 67media_77d= mean(Peso_77d); media_77d## [1] 845.9375\[H_O: \mu_1 = \mu_2\] \[H_A: \mu_1 \neq \mu_2\]
(prueba_3=t.test(rend$Peso_45d,rend$Peso_77d,alternative = "t", mu =0, paired=T,conf.level=0.95))##
## Paired t-test
##
## data: rend$Peso_45d and rend$Peso_77d
## t = -70.276, df = 15, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -802.5624 -755.3126
## sample estimates:
## mean of the differences
## -778.9375ifelse(prueba_3$p.value<0.05,"Rechazo Ho", " No Rechazo Ho")## [1] "Rechazo Ho"par(mfrow=c(1,2))
plot (rend$Peso_45d, xlab="dato",ylab= "peso del tuberculo (g)",main="Peso seco del tuberculo a los 45 dias", col="red",pch=16)
plot (rend$Peso_77d, xlab="dato",ylab= "peso del tuberculo (g)",main="Peso seco del tuberculo a los 77 dias", col="blue", pch=16)
cambio_rel =((media_45d-media_77d)*100)/media_77d; cambio_rel## [1] -92.07979cor(Peso_45d,Peso_77d,method = "pearson")## [1] 0.391328Ejercicio 3
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)
dataf=data.frame(Palma,Maiz); dataf## Palma Maiz
## 1 3 3
## 2 4 4
## 3 3 4
## 4 4 4
## 5 4 4
## 6 3 4
## 7 3 3
## 8 4 4
## 9 4 3
## 10 3 4
## 11 4 4
## 12 4 4
## 13 2 4
## 14 4 3
## 15 3 4
## 16 4 4
## 17 3 4
## 18 3 3
## 19 3 3
## 20 4 4
## 21 4 3par(mfrow=c(1,2))
median(Palma)## [1] 4median(Maiz)## [1] 4plot(Palma, col="Green",pch=16, xlab="dato",ylab="calidad de frito", main= "calidad de frito de las hojuelas de papa")
plot(Maiz, col="Brown",pch=12, xlab="dato",ylab="calidad de frito", main= "calidad de frito de las hojuelas de papa")
\[ H_0: Median_{palma} = Median_{Maiz} \\ H_a: Median_{palma} \neq Median_{Maiz} \]
Prueba_4=wilcox.test(dataf$Palma,dataf$Maiz, alternative="t", mu=0,paired =F, conf.level = 0.95)## Warning in wilcox.test.default(dataf$Palma, dataf$Maiz, alternative = "t", :
## cannot compute exact p-value with tiesifelse(Prueba_4$p.value<0.05, "Rechazo Ho", "Acepta Ho")## [1] "Acepta Ho"library(openxlsx)
excel=read.xlsx("C:/Users/faile/Desktop/aceite_diseño1.xlsx"); excel## L4C a4C b4C L12C a12C b12C
## 1 69.26 -1.31 28.68 62.20 0.81 37.31
## 2 68.15 -1.25 27.66 60.45 0.78 35.90
## 3 69.17 -1.49 28.02 63.12 0.55 36.36
## 4 68.88 -1.35 27.66 61.64 0.81 36.12
## 5 70.01 -1.32 27.66 61.25 0.77 36.45
## 6 70.15 -1.15 26.88 62.55 0.69 35.99
## 7 70.66 -1.25 26.25 64.12 0.59 36.14
## 8 68.68 -1.29 26.26 65.65 0.55 36.14
## 9 71.00 -1.42 28.15 66.87 0.42 35.55
## 10 72.18 -1.45 30.00 65.11 0.39 34.77
## 11 69.15 -1.29 28.24 66.14 0.41 32.32
## 12 70.00 -1.22 25.59 62.64 0.37 31.96
## 13 68.64 -1.19 24.69 61.97 0.35 30.17
## 14 68.12 -1.25 25.56 60.58 0.34 36.65
## 15 68.12 -1.25 26.26 60.68 0.34 37.15medianL4=median(excel$L4C); medianL4## [1] 69.17median_a4=median(excel$a4C);median_a4## [1] -1.29medianb4=median(excel$b4C);medianb4## [1] 27.66medianL12=median(excel$L12C);medianL12## [1] 62.55median_a12=median(excel$a12C);median_a12## [1] 0.55medianb12=median(excel$b12C);medianb12## [1] 36.12\[H_0:median_{L4C} =median{L12C} \\ H_0:median_{L4C} =median_{L12C}\]
wilcox_L=wilcox.test(excel$L4C,excel$L12C,alternative="t",mu=0, paired=T, conf.level=0.95)
ifelse(wilcox_L$p.value<0.05,"Rechaza Ho","Acepta Ho")## [1] "Rechaza Ho"\[H_0:median_{a4C} =median{a12C} \\ H_0:median_{a4C} =median_{a12C}\]
wilcox_a=wilcox.test(excel$a4C,excel$a12C,alternative="t",mu=0, paired=T, conf.level=0.95)## Warning in wilcox.test.default(excel$a4C, excel$a12C, alternative = "t", :
## cannot compute exact p-value with tiesifelse(wilcox_a$p.value<0.05,"Rechaza Ho","Acepta Ho")## [1] "Rechaza Ho"\[H_0:median_{b4C} =median_{b12C} \\ H_0:median_{b4C} =median_{b12C}\]
wilcox_b=wilcox.test(excel$b4C,excel$b12C,alternative="t",mu=0, paired=T, conf.level=0.95)
ifelse(wilcox_b$p.value<0.05,"Rechaza Ho","Acepta Ho")## [1] "Rechaza Ho"delta_e4=c(sqrt((excel$L4C)^2+(excel$a4C)^2+(excel$b4C)^2));delta_e4## [1] 74.97470 73.55991 74.64469 74.23848 75.28757 75.13241 75.38873 73.54042
## [9] 76.39004 78.17963 74.70532 74.54084 72.95520 72.76820 73.01702delta_e12=c(sqrt((excel$L12C)^2+(excel$a12C)^2+(excel$b12C)^2));delta_e12## [1] 72.53642 70.31089 72.84563 71.44788 71.27944 72.16827 73.60586 74.94214
## [9] 75.73358 73.81339 73.61556 70.32317 68.92483 70.80448 71.14985median(delta_e4)## [1] 74.64469median(delta_e12)## [1] 72.16827\[H_0: median_{\varDelta4C}=median_{\varDelta12C} \\ H_a: median_{\varDelta4C}\neq median_{\varDelta12C}\]
wilcox_delta_par= wilcox.test(delta_e4,delta_e12,alternative="t", mu=0, paired= T, conf.level=0.95)
ifelse(wilcox_delta_par$p.value<0.05,"Rechazo Ho", "Acepto Ho")## [1] "Rechazo Ho"library(scatterplot3d)scatterplot3d(excel[1:15,1:3],pch=16,color="Blue",main="Comportamiento del color a 4°C",col.grid = "black", type="h")
scatterplot3d(excel[1:15,4:6],pch=16,color = "red",main="Comportamiento del color a 12°C", col.grid="black",type="h")
# Ejercici 6: Comparar dos prevalencias
set.seed(1249)
# 0= sano, 1= enferma
xy=expand.grid(x=seq(1,12,1),y=seq(1,12,1))
lote1=round(runif(144,0,0.635),0)
color=ifelse(lote1==0,"green","red")
table(lote1)## lote1
## 0 1
## 119 25set.seed(1249)
# 0= sano, 1= enferma
xy_2=expand.grid(x=seq(1,12,1),y=seq(1,12,1))
lote2=round(runif(144,0,0.711),0)
color_2=ifelse(lote2==0,"green","red")
table(lote2)## lote2
## 0 1
## 105 39par(mfrow=c(1,2))
plot(xy,col=color,pch=16,main="Ditribucion de palmas en el lote 1")
plot(xy_2,col=color_2,pch=16,main="Ditribucion de palmas en el lote 2")
Calculo de las prevalencias \[Prevalencia=\frac{No.~casos~positivos}{Total ~de ~ casos} 100\]
table(lote1)## lote1
## 0 1
## 119 25table(lote2)## lote2
## 0 1
## 105 39p1=(25/144)*100;p1## [1] 17.36111p2=(39/144)*100;p2## [1] 27.08333\[Ho:prev_1=prev_2 \\ Ha: prev_1 \neq prev_2\]
(prueba_prev=prop.test(x=c(28,39),n=c(144,144)))# x=casos favorables, n= casos totales##
## 2-sample test for equality of proportions with continuity correction
##
## data: c(28, 39) out of c(144, 144)
## X-squared = 1.945, df = 1, p-value = 0.1631
## alternative hypothesis: two.sided
## 95 percent confidence interval:
## -0.18052771 0.02774993
## sample estimates:
## prop 1 prop 2
## 0.1944444 0.2708333ifelse(prueba_prev$p.value<0.05, "Rechazo Ho","Acepto Ho")## [1] "Acepto Ho"Ejercicio 7: tasas de incidencia
lote 1: 2 meses despues se registran 12 plantas enfermas adicionales
set.seed(1249)
# 0= sano, 1= enferma
xy1_nuevo=expand.grid(x=seq(1,12,1),y=seq(1,12,1))
lote1_nuevo=round(runif(144,0,0.69),0)
color_1nuevo=ifelse(lote1_nuevo==0,"green","red")
table(lote1_nuevo)## lote1_nuevo
## 0 1
## 107 37lote 2: 2 meses despues se registran 18 plantas enfermas adicionales
set.seed(1249)
# 0= sano, 1= enferma
xy_2nuevo=expand.grid(x=seq(1,12,1),y=seq(1,12,1))
lote2_nuevo=round(runif(144,0,0.885),0)
color_2_nuevo=ifelse(lote2_nuevo==0,"green","red")
table(lote2_nuevo)## lote2_nuevo
## 0 1
## 88 56par(mfrow=c(1,2))
plot(xy1_nuevo,col=color_1nuevo,pch=16,main="Ditribucion en el lote 1\n\ dos meses despues")
plot(xy_2nuevo,col=color_2_nuevo,pch=16,main="Ditribucion en el lote 2\n\ dos meses despues")
\[\lambda_1=\frac{conteo_1}{tiempo_1}\] \[\lambda_2=\frac{conteo_2}{tiempo_2}\]
conteo_lote1=12
conteo_lote2=18
tiempo_lote1=60 # se asumen dos meses de 30 dias cada uno
tiempo_lote2=60
#tasa de plantas enfermas por dia
lambda_lote1 = conteo_lote1/tiempo_lote1; lambda_lote1 ## [1] 0.2lambda_lote2 = conteo_lote2/tiempo_lote2; lambda_lote2## [1] 0.3#Tasas Poisson \[H_0:\lambda_1=\lambda_2\]
library(rateratio.test)
prueba_tasas =rateratio.test(c(12,18),c(tiempo_lote1,tiempo_lote2))# primer valor: conteos, segundo vector: magnitud tomada como referencia para las tasas
ifelse(prueba_tasas$p.value<0.05,"Rechazo Ho", "No rechazo Ho")## [1] "No rechazo Ho"# No hubo diferencias en los dos lotes, los dos genotipos estadisticamente presentan la misma tasa de incidencia de la enfermedad
Ejercicio 8: Prueba F- AOV-FSCA-B
# base de datos
library(readxl)
datos_ANOVA <- read_excel("C:/Users/faile/Desktop/Trabajos diseño/ANOVA taller_diseño.xlsx",
col_types = c("numeric", "numeric","numeric"))## Warning in read_fun(path = enc2native(normalizePath(path)), sheet_i = sheet, :
## Coercing text to numeric in A2 / R2C1: '7.1'## Warning in read_fun(path = enc2native(normalizePath(path)), sheet_i = sheet, :
## Coercing text to numeric in C2 / R2C3: '9.1'## Warning in read_fun(path = enc2native(normalizePath(path)), sheet_i = sheet, :
## Coercing text to numeric in A13 / R13C1: '7.4'## Warning in read_fun(path = enc2native(normalizePath(path)), sheet_i = sheet, :
## Coercing text to numeric in C13 / R13C3: '8.7'View(datos_ANOVA)# factores
metodo= c(rep("Bray",22),rep("Olsen",22),rep("Mehlich",22)); metodo## [1] "Bray" "Bray" "Bray" "Bray" "Bray" "Bray" "Bray"
## [8] "Bray" "Bray" "Bray" "Bray" "Bray" "Bray" "Bray"
## [15] "Bray" "Bray" "Bray" "Bray" "Bray" "Bray" "Bray"
## [22] "Bray" "Olsen" "Olsen" "Olsen" "Olsen" "Olsen" "Olsen"
## [29] "Olsen" "Olsen" "Olsen" "Olsen" "Olsen" "Olsen" "Olsen"
## [36] "Olsen" "Olsen" "Olsen" "Olsen" "Olsen" "Olsen" "Olsen"
## [43] "Olsen" "Olsen" "Mehlich" "Mehlich" "Mehlich" "Mehlich" "Mehlich"
## [50] "Mehlich" "Mehlich" "Mehlich" "Mehlich" "Mehlich" "Mehlich" "Mehlich"
## [57] "Mehlich" "Mehlich" "Mehlich" "Mehlich" "Mehlich" "Mehlich" "Mehlich"
## [64] "Mehlich" "Mehlich" "Mehlich"fosforo=c(datos_ANOVA$Bray,datos_ANOVA$Olsen,datos_ANOVA$Mehlich);fosforo## [1] 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.0 6.9 6.7 7.1
## [20] 6.6 6.6 6.9 6.0 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
## [39] 6.6 6.3 5.8 6.5 6.2 6.4 9.1 7.1 7.8 7.3 7.6 7.8 7.4 7.3 7.3 7.1 8.0 8.7 7.6
## [58] 7.8 7.8 7.2 7.2 7.9 7.3 7.5 7.1 5.1medidas descriptivas
# medias
medias_anova=apply(datos_ANOVA,2,mean);medias_anova## Bray Olsen Mehlich
## 6.813636 6.409091 7.500000# desviaciones
desviaciones_anova=apply(datos_ANOVA,2,sd);desviaciones_anova## Bray Olsen Mehlich
## 0.2030610 0.2408499 0.7361418# coef.variacion
cv_anova=(desviaciones_anova/medias_anova)*100;cv_anova## Bray Olsen Mehlich
## 2.980215 3.757941 9.815224# boxplot
boxplot(datos_ANOVA,datos_ANOVA,col=c(rep("lightblue",1),rep("lightgreen",1),rep("pink",1)),ylab="mg Fosforo / Kg de suelo",main=" Resultados de muestras de fosforo en\n\ suelo usando tres metodos distintos",cex.main=1)
points(c(1:3),medias_anova,col="red",pch=16)
# Analisis de varianza \[Ho: \tau_{Bray} = \tau_{Olsen} = \tau_{Mehlich}\]
ANOVA=aov(fosforo~metodo,datos_ANOVA)
resumen=summary(ANOVA); resumen## Df Sum Sq Mean Sq F value Pr(>F)
## metodo 2 13.38 6.691 31.31 3.62e-10 ***
## Residuals 63 13.46 0.214
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1p.valor=resumen[[1]][1,5];p.valor # accede a los objetos de la lista, y toma en especifico el objeto de la fila 1 y columna 5## [1] 3.624612e-10ifelse(p.valor<0.05,"Rechazo Ho","Acepto Ho")## [1] "Rechazo Ho"Estadisticamente el metodo utilizado tiene influencia en la extraccion de fosforo en las muestras de suelo
Comparacion de medias (prueba de Tukey)
TukeyHSD(ANOVA,'metodo')## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = fosforo ~ metodo, data = datos_ANOVA)
##
## $metodo
## diff lwr upr p adj
## Mehlich-Bray 0.6863636 0.3517897 1.02093754 0.0000190
## Olsen-Bray -0.4045455 -0.7391194 -0.06997155 0.0139314
## Olsen-Mehlich -1.0909091 -1.4254830 -0.75633518 0.0000000los tratamientos son muy diferentes el uno con el otro # Normalidad de los residuales \[Ho: \epsilon_{ij}\sim N(0,\sigma^2_e)\]
residuales=residuals(ANOVA)
shapiro.test(residuales)##
## Shapiro-Wilk normality test
##
## data: residuales
## W = 0.75936, p-value = 4.772e-09*los residuales no tienen un comportamiento ajustado a la normal
prueba homocedasticidad (igualdad de varianzas)
\[Ho: \sigma^2_{Bray} = \sigma^2_{Olsen} = \sigma^2_{Mehlich}\]
homocedasticidad=bartlett.test(fosforo~metodo)
ifelse(homocedasticidad$p.value<0.05,"Rechaza Ho","Acepta Ho")## [1] "Rechaza Ho"