graphics.off()
cat("\014")
require(table1)
## Cargando paquete requerido: table1
## Warning: package 'table1' was built under R version 4.5.3
## 
## Adjuntando el paquete: 'table1'
## The following objects are masked from 'package:base':
## 
##     units, units<-
require(ggplot2)
## Cargando paquete requerido: ggplot2
require(agricolae)
## Cargando paquete requerido: agricolae
## Warning: package 'agricolae' was built under R version 4.5.3
require(car)
## Cargando paquete requerido: car
## Cargando paquete requerido: carData
load("Salinidad.RData")
str(Salinidad)
## 'data.frame':    45 obs. of  5 variables:
##  $ Biomasa  : num  765 954 828 755 896 ...
##  $ pH       : num  5 4.7 4.2 4.4 5.55 5.5 4.25 4.45 4.75 4.6 ...
##  $ Salinidad: int  33 35 32 30 33 33 36 30 38 30 ...
##  $ Zinc     : num  16.5 14 15.3 17.3 22.3 ...
##  $ Potasio  : num  1442 1299 1154 1045 522 ...
#Análisis exploratorio univariado para cada característica
sapply(Salinidad, mean)
##     Biomasa          pH   Salinidad        Zinc     Potasio 
## 1082.172644    4.608889   30.266667   17.830796  797.377778
sapply(Salinidad, sd)
##    Biomasa         pH  Salinidad       Zinc    Potasio 
## 546.287433   1.254731   3.719726   8.274169 297.576022
cv = sapply(Salinidad, sd) / sapply(Salinidad, mean)*100
round(cv,1)
##   Biomasa        pH Salinidad      Zinc   Potasio 
##      50.5      27.2      12.3      46.4      37.3
table1(~Biomasa + pH + Salinidad + Zinc + Potasio, data=Salinidad)
Overall
(N=45)
Biomasa
Mean (SD) 1080 (546)
Median [Min, Max] 992 [370, 2340]
pH
Mean (SD) 4.61 (1.25)
Median [Min, Max] 4.45 [3.20, 7.45]
Salinidad
Mean (SD) 30.3 (3.72)
Median [Min, Max] 30.0 [24.0, 38.0]
Zinc
Mean (SD) 17.8 (8.27)
Median [Min, Max] 19.2 [0.211, 31.3]
Potasio
Mean (SD) 797 (298)
Median [Min, Max] 773 [351, 1440] 
#------------BIOMASA---------
#La biomasa presentó una media de 1080 g y una desviación estándar de 546 g,
#con valores entre 370 y 2340 g. El CV fue de 50.5% indicando una variabilidad
#relativamente alta entre las muestras. Mediana fue de 992 g.
#----------pH-----------------
#El pH presentó una media de 4.61, desviación estándar de 1.25 y valores entre 3.2 y 7.45. Su
#coeficiente de variación fue de 27.2%. Mediana de 4.45.
#----------------SALINIDAD------------------
#La salinidad presentó una media de 30.3, desviación estándar de 3.72 y un rango entre 24 y 38.
#El CV fue de 12.3%, siendo la variable con menor variabilidad relativa. Mediana de 30.
#------------------ZINC----------------
#Presentó una media de 17.8, desviación estándar de 8.27 y valores entre 0.2011 y 31.3. Su CV fue de 46.4%,
#indicando una alta variabilidad relativa. Mediana de 19.2.
#------------POTASIO-----------------------
#El potasio presentó una media de 797, desviación estándar de 298 y valores entre 351 y 1440. El CV
#fue de 37.3%, por lo que presentó una variabilidad relativa intermedia-alta. Mediana de 773.
par(mfrow = c(2,3))
hist(Salinidad$Biomasa, main="Biomasa", xlab="gr")
#Media mayor a la mediana. Distribución asimétrica a la derecha (positiva). Una moda principal en 500-1500 y un grupito aparte en 2000-2500.
hist(Salinidad$pH,        main="pH",        xlab="pH")
#Media ligeramente mayor a la mediana, aunque la presencia de valores elevados entre aproximadamente
# 7.1 y 7.45 incrementa la dispersión y genera cierta asimetría.Asimetría positiva. Varias modas: 
# una en 3-3.5, otra en 4.5-5 y un grupo en 7+.
hist(Salinidad$Salinidad, main="Salinidad", xlab="Salinidad")
#Media prácticamente igual a la mediana. Casi simétrica. varias modas: una en 24-26, otra en 28-30 y en 34-36 (bimodal).
hist(Salinidad$Zinc,      main="Zinc",      xlab="Zinc")
#Media menor a la mediana. Asimetría negativa. Una moda en 15-25, con un grupito aparte cerca de 0.
hist(Salinidad$Potasio,   main="Potasio",   xlab="Potasio")
#Media un poco mayor a la mediana. Asimetría positiva (leve). Varias modas: 400-600, 800-1000 y 1200-1400.
par(mfrow = c(2,3))

boxplot(Salinidad$Biomasa,main="Biomasa")
#La mediana está casi en el centro de la caja, pero el bigote de arriba es mucho más largo que el de abajo.
#Indica asimetría positiva.
#No hay valores atípicos.Las muestras más altas están en el extremo del bigote.
boxplot(Salinidad$pH,main="pH")
#La caja va de 3.45 a 5.35 y la mdiana (4.45) está ubicada en medio. El bigote inferior es corto y el superior largo.
#Asimetría positiva.
#No hay valores atípicos.
boxplot(Salinidad$Salinidad, main="Salinidad")
#La mediana (30) queda justo en el centro de la caja y los bigotes son parecidos, el largo siendo un poco más largo.
#Es casi simétrica.
#Sin atípicos.
boxplot(Salinidad$Zinc,main="Zinc")
#Mediana más cerca de Q3 que de Q1, bigote de abajo termina en 9.4. Sumado al punto suelto, indica asimetría negativa.
#Valor atípico que representa a cinco muestras (0.21 a 0.37). Forman un grupo con características propias y no errores de medición.
boxplot(Salinidad$Potasio, main="Potasio")
#Bigote superior mucho más largo que el inferior. Asimetría positiva leve. Dentro de la caja, mediana un poco más cerca de Q3.
#así que las señales no son del todo claras.
#Sin atípicos.
par(mfrow = c(1,1))

#Análisis exploratorio bivariado
R2 = cor(Salinidad[, c("Biomasa", "pH", "Salinidad", "Zinc", "Potasio")])^2
round(R2, 3)
##           Biomasa    pH Salinidad  Zinc Potasio
## Biomasa     1.000 0.861     0.004 0.611   0.005
## pH          0.861 1.000     0.002 0.519   0.001
## Salinidad   0.004 0.002     1.000 0.182   0.000
## Zinc        0.611 0.519     0.182 1.000   0.006
## Potasio     0.005 0.001     0.000 0.006   1.000
#pH: explica casi toda la variación de la biomasa. Queda un 13.9% sin explicar, que es la dispersión de los puntos alrededor de la recta.
#Zinc: también explica mucho, pero unos 25 puntos porcentuales menos que el pH.
#Salinidad y potasio: explican menos del 1%. Con 45 muestras, un R^2 tan bajo no es significativo. No sirve para predecir la biomasa. 
ggplot(Salinidad, aes(x=pH, y=Biomasa)) + geom_point() + geom_smooth(method="lm")
## `geom_smooth()` using formula = 'y ~ x'

#Los puntos siguen la recta de cerca y el intervalo de confianza del 95% es angosta, es decir, recta bien estimada de poca incertidumbre.
#A mayor pH, mayor biomasa.
#El pH explica cerca del 86% de la biomasa (R^200.86), este dato se confirma más abajo.
ggplot(Salinidad, aes(x=Zinc, y=Biomasa)) + geom_point() + geom_smooth(method="lm")
## `geom_smooth()` using formula = 'y ~ x'

#La relación es negativa: a más Zinc, menos biomasa. Su R^2=0.61 es fuerte, pero más dispersa que la del pH.
#Entre Zinc 10 y 31 los puntos bajan con más dispersión, y uno con Zinc ^20 y biomasa ~1890 destaca por encima.
#Parte de esta relación puede ser indirecta: pH y Zinc están correlacionados entre sí (R^2=-0.72), y son las mismas muestras
#en el extresmo. No se puede decir que el Zinc reduzca la biomasa. Solo que ambas varían juntas.
ggplot(Salinidad, aes(x=Salinidad, y=Biomasa)) + geom_point() + geom_smooth(method="lm")
## `geom_smooth()` using formula = 'y ~ x'

#La recta es casi horizontal y el intervalo de confianza es muy ancho.No se distingue el cero.
#Los puntos forman columnas verticales porque la salinidad solo toma valores enteros. En cada columna hay biomasas muy distintas: con
#salinidad 30, van de unos 370 a más de 2300 g.
#La salinidad no explica la biomasa en este conjunto de datos.
ggplot(Salinidad, aes(x=Potasio, y=Biomasa)) + geom_point() + geom_smooth(method="lm")
## `geom_smooth()` using formula = 'y ~ x'

#Tampoco hay tendencia: recta casi plana y banda ancha que incluye la horizontal.
#Se ven dos zonas: 5 puntos arriba a la izquierda (Potasio ~450-570, biomasa mayor a 2100) y muchos abajo con 
#biomasa baja y Potasio parecido. Con Potasio ~500, la biomasa puede ser de 370 a 2337 g, así que el Potasio no 
#sirve para predecirla.
#Categorizar, ANOVA y LSD

if ("pH_cat" %in% names(Salinidad)) Salinidad$pH_cat = NULL   

cortes = quantile(Salinidad$pH, probs = c(0, 1/3, 2/3, 1))
cortes
##        0% 33.33333% 66.66667%      100% 
##  3.200000  3.883333  4.900000  7.450000
cortes = quantile(Salinidad$pH, probs = c(0, 1/3, 2/3, 1))
cortes
##        0% 33.33333% 66.66667%      100% 
##  3.200000  3.883333  4.900000  7.450000
pH_cat = cut(Salinidad$pH, breaks = cortes, include.lowest = TRUE,
             labels = c("Bajo", "Medio", "Alto"))
table(pH_cat)
## pH_cat
##  Bajo Medio  Alto 
##    15    15    15
#El pH se dividió en terciles porque da grupos de tamaño exactamente igual (15 cada uno), sin imponer cortes arbitrarios:
#Bajo(3.20 a 3.88), Medio(3.88 a 4.90) y Alto(4.90 a 7.45).
table1(~Biomasa|pH_cat, data = data.frame(Biomasa = Salinidad$Biomasa, pH_cat = pH_cat))
Bajo
(N=15)		 Medio
(N=15)		 Alto
(N=15)		 Overall
(N=45)
Biomasa
Mean (SD) 593 (165) 1050 (245) 1610 (548) 1080 (546)
Median [Min, Max] 546 [370, 978] 1040 [568, 1490] 1420 [765, 2340] 992 [370, 2340] 
#La biomasa media aumenta escalonadamente con el nivel de pH.
#En Bajo y Medio la media y la mediana casi coinciden (distribuciones simétricas); en Alto la media (1610) supera bastante 
#a la mediana (1420), señal de asimetría positiva, coherente con las 5 muestras extremas (biomasa > 2100) que tiene este grupo.
ggplot(Salinidad, aes(x=pH_cat, y=Biomasa, fill=pH_cat)) + geom_boxplot()

#En el boxplot, las tres cajas casi no se tralapan y suben de izquierda a derecha. Bajo tiene un atípico leve (~978 g). La caja 
#de Alto es mucho más ancha que las otras dos: va de ~1200 a ~2190, contra ~480-660 en Bajo y ~890-1200 en Medio. Esa diferencia
#de tamaño es la primera pista visual de que las varianzas no son iguales.
anova_sal = aov(Salinidad$Biomasa ~ pH_cat)
#H0: las medias de biomasa son iguales en los tres niveles de pH. H1: al menos un nivel difiere.
#Como p<0.05, se rechaza H0: el nivel de pH genera diferencias significativas en la biomasa.
summary(anova_sal)
##             Df  Sum Sq Mean Sq F value   Pr(>F)    
## pH_cat       2 7712683 3856342   29.89 8.45e-09 ***
## Residuals   42 5418235  129006                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
shapiro.test(residuals(anova_sal))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(anova_sal)
## W = 0.97316, p-value = 0.3749
#H0: los residuales siguen una distribución normal, Como p>0.05, no s rechaza H0: el supuesto se cumple. 
qqnorm(residuals(anova_sal)); qqline(residuals(anova_sal))

#El Q-Qplot lo confirma: la mayoría de los puntos, sobre todo en el centro, sigue de cerca la recta. hay una desviación leve en ambas colas 
#algunos puntos por debajo de la recta a la izquierda y por encima a la derecha), lo que indica colas un poco más pesadas que una normal perfecta, pero no
#lo suficiente como para que Shapiro lo detecte como anormal. Es un resultado coherente y no contradictorio.
bartlett.test(Salinidad$Biomasa ~ pH_cat)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  Salinidad$Biomasa by pH_cat
## Bartlett's K-squared = 20.084, df = 2, p-value = 4.353e-05
leveneTest(Salinidad$Biomasa ~ pH_cat)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value   Pr(>F)    
## group  2  8.6753 0.000702 ***
##       42                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#H0: las tres varianzas son iguales. Ambas pruebas dan p<0.05, así que se rechaza H0 en las dos: el supuesto no se cumple. 
#Que las dos pruebas coincidan le da más fuerza a la conclusión: no es un artefacto de una sola prueba.
#La causa se ve en las DE: 165(Bajo), 245(Medio) y 548(Alto). El grupo Alto es 3.3 veces más variable que bajo, porque mezcla pH "normales" (5-5.6) con el grupo 
#aparte de pH mayor a 7 (biomasa 2160.2337).
oneway.test(Salinidad$Biomasa ~ pH_cat)
## 
##  One-way analysis of means (not assuming equal variances)
## 
## data:  Salinidad$Biomasa and pH_cat
## F = 34.727, num df = 2.000, denom df = 24.736, p-value = 6.58e-08
#El resultado confirma al ANOVA clásico que el pH afecta la biomasa.
LSD.test(anova_sal, "pH_cat", console = TRUE)
## 
## Study: anova_sal ~ "pH_cat"
## 
## LSD t Test for Salinidad$Biomasa 
## 
## Mean Square Error:  129005.6 
## 
## pH_cat,  means and individual ( 95 %) CI
## 
##       Salinidad.Biomasa      std  r      se       LCL       UCL     Min
## Alto          1605.3853 547.6014 15 92.7382 1418.2321 1792.5386 765.280
## Bajo           593.0233 164.8299 15 92.7382  405.8701  780.1766 369.823
## Medio         1048.1093 244.9093 15 92.7382  860.9560 1235.2625 568.455
##            Max       Q25      Q50      Q75
## Alto  2337.326 1200.2575 1422.836 2192.560
## Bajo   977.515  481.3520  545.538  659.713
## Medio 1491.276  890.8515 1039.637 1198.396
## 
## Alpha: 0.05 ; DF Error: 42
## Critical Value of t: 2.018082 
## 
## least Significant Difference: 264.6747 
## 
## Treatments with the same letter are not significantly different.
## 
##       Salinidad$Biomasa groups
## Alto          1605.3853      a
## Medio         1048.1093      b
## Bajo           593.0233      c
#Letras distintas para los tres grupos dignifica que los tres niveles de pH difieren significativamente entre sí.

R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

summary(cars)
##      speed           dist       
##  Min.   : 4.0   Min.   :  2.00  
##  1st Qu.:12.0   1st Qu.: 26.00  
##  Median :15.0   Median : 36.00  
##  Mean   :15.4   Mean   : 42.98  
##  3rd Qu.:19.0   3rd Qu.: 56.00  
##  Max.   :25.0   Max.   :120.00

Including Plots

You can also embed plots, for example:

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.