Este documento contiene los códigos base que necesitas para realizar el análisis de los datos de la investigación sobre la variabilidad poblacional (dentro y entre poblaciones de semillas). Debes completer aquellos cálculos y/o análisis que se indiquen y redactar tu propia Introducción, Metadata, Descripción de Resultados, Discusión, Conclusiones y Bibliografía.
Escribir uno o dos párrafos de información de trasfondo sobre la
especie (su reproducción) y variabilidad poblacional en general y sus
correlativos (posibles causas):
* No más de 500 palabras
* Con al menos dos citas (incluya en el texto: autor, año) de
publicaciones científicas.
Formular preguntas e hipótesis sobre los siguientes aspectos:
Variación en las medidas (estadísticas descriptivas) de las variables medidas (longitud, ancho máximo, masa), dentro de cada población (árbol).
Distribución de los valores de las variables para cada árbol y en conjunto.
Diferencias entre las poblaciones (árboles) en la longitud ancho máximo y masa de las semillas.
Diferencias entre las poblaciones en el número de semillas por fruto.
Relación entre el número de semillas por fruto y la longitud y masa de las semillas.
Relación entre la masa y la longitud (o ancho máximo) de las semillas.
Las hipótesis deben contener una premisa (dado que…) y la predicción (espero que…).
En este chunk debes incluir la activación de los paquetes que necesita. Para activar los paquetes, estos deben estar previamente instalados.
library(ggplot2)
library(dplyr)
library(EnvStats)
library(gt)
library(agricolae)
library(Hmisc)
Los datos completos deben estar en la estructura indicada en el laboratorio y en un archivo CSV.
semillas <- read.csv("variacion poblacional 2023.csv")
head(semillas)
## grupo arbol fruto semilla long.mm amax.mm masa.g n.sem.fruto
## 1 1 1 1 1 11.3 8.6 0.3290 4
## 2 1 1 1 2 11.4 6.0 0.2355 4
## 3 1 1 1 3 11.2 8.7 0.3074 4
## 4 1 1 1 4 12.2 6.5 0.2988 4
## 5 1 1 2 1 11.2 7.4 0.3233 6
## 6 1 1 2 2 12.2 10.3 0.3483 6
Aquí debes escribir la metadata (ver detalles en
MOODLE) del proyecto, la cual debe incluir:
- Título del proyecto
Fecha de la toma de datos
Participantes
Lugar
Coordenadas
Hora de inicio y final de las mediciones
Organismos estudiados (nombre científico)
Diseño de muestreo y toma de datos
Métodos de mediciones
Lista de variables o atributos, con la siguiente información:
Después de cada sección de resultados debe escribir un párrafo en el que describa los resultados más notables y relevantes.
Con las estadísticas descriptivas podemos empezar a responder las preguntas sobre las variaciones en las variables medidas, dentro y entre poblaciones (árboles).
En esta sección calculamos media (promedio), mediana, desviación y error estándar, sesgo y curtosis (debes buscar información sobre los estadísticos que no conoces).
Estos estadísticos nos darán información sobre la tendencia central y la dispersión de los datos en cada población (árbol). Los cálculos debes repetirlos para ancho máximo y masa. Estos resultados también pueden ayudar a detectar errores en los datos.
#función para definir estadísticos a calcular
mystats <- function(x){
m <- mean(x)
md <- median(x)
n <- length(x)
s <- sd(x)
se <- sd(x)/sqrt(length(x))
skew <- sum((x-m)^3/s^3)/n
kurt <- sum((x-m)^4/s^4)/n - 3
return(c(n=n, mean=m, median=md, stdev=s, serror=se, skew=skew, kurtosis=kurt))
}
# función sapply con split para aplicar la función a datos de cada árbol
tabla.estad <- sapply(split(semillas$long.mm, semillas$arbol), mystats)
colnames(tabla.estad) <- c("Arbol1", "Arbol2", "Arbol3")
rownames(tabla.estad) <- c("n", "Media", "Mediana", "Desv. Est.", "Error Est.", "Sesgo", "Curtosis")
tabla.estad
## Arbol1 Arbol2 Arbol3
## n 47.00000000 58.0000000 81.00000000
## Media 11.80851064 10.0086207 9.13481481
## Mediana 11.80000000 9.9000000 9.10000000
## Desv. Est. 0.77201875 0.8316723 0.49745379
## Error Est. 0.11261051 0.1092039 0.05527264
## Sesgo -0.05021799 0.7969947 0.06599508
## Curtosis -0.76310913 0.2042062 -0.30422115
Utilizando el paquete gt se pueden crear tablas más elegantes. Las tablas deben llevar una leyenda (no es título) en la parte superior.
table1 <- gt(
as.data.frame(tabla.estad),
rownames_to_stub = TRUE
) %>%
tab_header("TABLA 1. Estadísticas descriptivas de la longitud de semillas en cada población"
) %>%
tab_stubhead(label = "Estadístico"
) %>%
cols_label(
Arbol1 = md("Árbol 1"),
Arbol2 = md("Árbol 2"),
Arbol3 = md("Arbol 3")
) %>%
fmt_number(
columns = c(Arbol1,Arbol2,Arbol3),
decimals = 2)
table1
TABLA 1. Estadísticas descriptivas de la longitud de semillas en cada población | |||
Estadístico | Árbol 1 | Árbol 2 | Arbol 3 |
---|---|---|---|
n | 47.00 | 58.00 | 81.00 |
Media | 11.81 | 10.01 | 9.13 |
Mediana | 11.80 | 9.90 | 9.10 |
Desv. Est. | 0.77 | 0.83 | 0.50 |
Error Est. | 0.11 | 0.11 | 0.06 |
Sesgo | −0.05 | 0.80 | 0.07 |
Curtosis | −0.76 | 0.20 | −0.30 |
Una manera de visualizar los datos para ver la tendencia central y dispersión, en cada población, es mediante el uso de histogramas. También podemos incorporar a la gráfica la curva de distribución normal, basada en los estadísticos media y desviación estándar. Repetir para las otras variables.
# histograma y curva normal para cada arbol; variable longitud
# árbol 1
arbol1 <- semillas %>%
filter(arbol == "1")
msa1 <- semillas %>%
filter(arbol == "1") %>%
summarise(means = mean(long.mm), sd=sd(long.mm))
ggplot(arbol1, aes(x=long.mm)) +
geom_histogram(aes(y=..density..), position="identity", binwidth=.25, color="blue", fill="gray") +
stat_function(fun = dnorm, color="blue", size=1, args=list(mean=msa1$means, sd=msa1$sd)) +
labs(x="Longitud semilla, mm", y="densidad") +
annotate("text", x = 11, y = 1, label = "Arbol 1")
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
# árbol 2
arbol2 <- semillas %>%
filter(arbol == "2")
msa2 <- semillas %>%
filter(arbol == "2") %>%
summarise(means = mean(long.mm), sd=sd(long.mm))
ggplot(arbol2, aes(x=long.mm)) +
geom_histogram(aes(y=..density..), position="identity", binwidth=.25, color="orange", fill="gray") +
stat_function(fun = dnorm, color="orange", size=1, args=list(mean=msa2$means, sd=msa2$sd)) +
labs(x="Longitud semilla, mm", y="densidad") +
annotate("text", x = 11, y = 1, label = "Arbol 2")
# arbol 3
arbol3 <- semillas %>%
filter(arbol == "3")
msa3 <- semillas %>%
filter(arbol == "3") %>%
summarise(means = mean(long.mm), sd=sd(long.mm))
ggplot(arbol3, aes(x=long.mm)) +
geom_histogram(aes(y=..density..), position="identity", binwidth=.25, color="green", fill="gray") +
stat_function(fun = dnorm, color="green", size=1, args=list(mean=msa3$means, sd=msa3$sd)) +
labs(x="Longitud semilla, mm", y="densidad") +
annotate("text", x = 11, y = 1, label = "Arbol 3")
Escribir leyenda debajo de cada gráfica: FIGURA 1. Histograma de la longitud de las semillas…
Mediante esta gráfica podemos comparar las tres poblaciones en una misma gráfica. Repetir para las otras variables.
#las tres curvas normales juntas
ggplot(semillas, aes(x=long.mm)) +
stat_function(fun = dnorm, args=list(mean=msa1$means, sd=msa1$sd), aes(colour = "Arbol 1"), size=1.5) +
stat_function(fun = dnorm, args=list(mean=msa2$means, sd=msa2$sd), aes(colour = "Arbol 2"), size=1.5) +
stat_function(fun = dnorm, args=list(mean=msa3$means, sd=msa3$sd), aes(colour = "Arbol 3"), size=1.5) +
scale_colour_manual("T. populnea", values = c("blue","orange","green")) +
labs(x="Longitud semilla, mm", y="densidad")
Repetir para las otras variables tomando en cuenta lo siguiente:
long <- ggplot(semillas, aes(x = arbol, y=long.mm, group=arbol)) + geom_boxplot(fill=c("blue", "orange", "green")) +
stat_summary(aes(group=arbol), fun=mean, colour="darkred", geom="point", shape=18, size=3) +
labs(x = "Arbol", y = "Longitud semilla, mm")
long
Repetir para las otras variables.
Las siguientes pruebas estadísticas se usan para determinar si nuestros datos cumplen con los supuestos para realizar pruebas paramétricas.
Debes repetirla para cada una de las variables.
#Q-Q plot
EnvStats::qqPlot(semillas$long.mm, add.line = TRUE, points.col = "blue", line.col = "red")
#prueba Shapiro-Wilk
shapiro.test(semillas$long.mm)
##
## Shapiro-Wilk normality test
##
## data: semillas$long.mm
## W = 0.92253, p-value = 2.304e-08
Esta es la prueba mas importante para decidir si podemos usar métodos paramétricos.
bartlett.test(long.mm ~ arbol, data=semillas)
##
## Bartlett test of homogeneity of variances
##
## data: long.mm by arbol
## Bartlett's K-squared = 19.551, df = 2, p-value = 5.683e-05
La prueba de ANOVA (análisis de la varianza) nos permite probar si hay un efecto significativo en la variable, debido a los factores incluidos.
anova_l <- aov(long.mm ~ as.factor(arbol) * as.factor(fruto), data=semillas)
summary(anova_l)
## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(arbol) 2 213.08 106.54 253.472 < 2e-16 ***
## as.factor(fruto) 2 2.95 1.48 3.509 0.032017 *
## as.factor(arbol):as.factor(fruto) 4 9.29 2.32 5.526 0.000325 ***
## Residuals 177 74.40 0.42
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Con esta prueba podemos discriminar cuáles son las parejas de árboles que son diferentes entre si.
anova_l_arbol <- aov(long.mm ~ as.factor(arbol), data=semillas)
summary(anova_l_arbol)
## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(arbol) 2 213.08 106.54 225 <2e-16 ***
## Residuals 183 86.64 0.47
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc <- TukeyHSD(anova_l_arbol)
posthoc
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = long.mm ~ as.factor(arbol), data = semillas)
##
## $`as.factor(arbol)`
## diff lwr upr p adj
## 2-1 -1.7998899 -2.118982 -1.4807978 0
## 3-1 -2.6736958 -2.971820 -2.3755712 0
## 3-2 -0.8738059 -1.153469 -0.5941427 0
par(cex = 1)
plot(posthoc)
Cuando se violan los supuestos del ANOVA, una alternativa es la prueba de Kruskal-Wallis, que no depende de la varianza.
kruskal.test(long.mm ~ arbol, data = semillas)
##
## Kruskal-Wallis rank sum test
##
## data: long.mm by arbol
## Kruskal-Wallis chi-squared = 119.86, df = 2, p-value < 2.2e-16
Para probar si existe una relación entre dos variables, podemos usar el análisis de correlación, el cual no implica que haya una relación causa-efecto entre las variables, solo una asociación. Si establecemos una hipótesis de causa-efecto, entonces podemos utilizar el análisis de regresión para probar esta relación (modelo).
Vamos a realizar un análisis de correlación para las variables, en conjunto y para cada árbol.
correla <- semillas |>
select(long.mm,amax.mm,masa.g)
corr_arboles <- rcorr(as.matrix(correla))
round(corr_arboles$r, digits = 4)
## long.mm amax.mm masa.g
## long.mm 1.0000 0.1307 0.6605
## amax.mm 0.1307 1.0000 0.4349
## masa.g 0.6605 0.4349 1.0000
round(corr_arboles$P, digits = 4)
## long.mm amax.mm masa.g
## long.mm NA 0.0754 0
## amax.mm 0.0754 NA 0
## masa.g 0.0000 0.0000 NA
correlaA1 <- semillas |>
filter(arbol == 1) |>
select(long.mm,amax.mm,masa.g)
corr_arbol1 <- rcorr(as.matrix(correlaA1))
round(corr_arbol1$r, digits = 4)
## long.mm amax.mm masa.g
## long.mm 1.0000 -0.1925 0.2664
## amax.mm -0.1925 1.0000 0.2394
## masa.g 0.2664 0.2394 1.0000
round(corr_arbol1$P, digits = 4)
## long.mm amax.mm masa.g
## long.mm NA 0.1949 0.0703
## amax.mm 0.1949 NA 0.1051
## masa.g 0.0703 0.1051 NA
correlaA2 <- semillas |>
filter(arbol == 2) |>
select(long.mm,amax.mm,masa.g)
corr_arbol2 <- rcorr(as.matrix(correlaA2))
round(corr_arbol2$r, digits = 4)
## long.mm amax.mm masa.g
## long.mm 1.0000 -0.2786 0.1456
## amax.mm -0.2786 1.0000 0.4257
## masa.g 0.1456 0.4257 1.0000
round(corr_arbol2$P, digits = 4)
## long.mm amax.mm masa.g
## long.mm NA 0.0342 0.2754
## amax.mm 0.0342 NA 0.0009
## masa.g 0.2754 0.0009 NA
correlaA3 <- semillas |>
filter(arbol == 3) |>
select(long.mm,amax.mm,masa.g)
corr_arbol3 <- rcorr(as.matrix(correlaA3))
round(corr_arbol3$r, digits = 4)
## long.mm amax.mm masa.g
## long.mm 1.0000 0.0684 0.1882
## amax.mm 0.0684 1.0000 0.5566
## masa.g 0.1882 0.5566 1.0000
round(corr_arbol3$P, digits = 4)
## long.mm amax.mm masa.g
## long.mm NA 0.5443 0.0925
## amax.mm 0.5443 NA 0.0000
## masa.g 0.0925 0.0000 NA
En el análisis de regresión planteamos un modelo lineal y probamos si este satisface una posible relación entre las dos variables, usando el número de semillas por fruto como la variable independiente y la masa de semillas individuales como la variable respuesta.
ggplot(semillas, aes(x = n.sem.fruto, y = masa.g)) +
geom_point() +
geom_smooth(method = "lm", se = TRUE)
## `geom_smooth()` using formula = 'y ~ x'
regres <- lm(masa.g ~ n.sem.fruto, data = semillas)
summary(regres)
##
## Call:
## lm(formula = masa.g ~ n.sem.fruto, data = semillas)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.14842 -0.03260 -0.01232 0.03735 0.16658
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.31252 0.01091 28.644 < 2e-16 ***
## n.sem.fruto -0.01070 0.00125 -8.562 4.33e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05089 on 184 degrees of freedom
## Multiple R-squared: 0.2849, Adjusted R-squared: 0.281
## F-statistic: 73.31 on 1 and 184 DF, p-value: 4.335e-15
En esta sección debe explicar sus resultados, de acuerdo a sus preguntas e hipótesis, y comparando con la literatura (al menos dos referencias de literatura científica primaria).
En esta sección debe seleccionar de los resultados y discusión aquello que constituye lo más relevante y general sobre la investigación. Debe incluir propuestas de nuevas investigaciones e hipótesis sobre el tema. También puede incluir reflexiones personales sobre la investigación y sugerencias al profesor.
Las referencias deben estar en un formato completo, usado en alguna publicación sobre ecología (Biotrópica, por ejemplo).