Bienvenid@s a este módulo de introducción a la sintaxis en estadística en R. Este módulo pretende simplemente extender los contenidos para que estéis familiarizados con algunos de los procedimientos de uso de la estadística en R. Pero no pretende ser un curso de estadística en sí mismo ya que proporcionar aquí contenidos teóricos sobre muestreos, estimadores, teoría de la probabilidad, etc. excede la duración del curso. Pero sí creemos que proporciona una gran utilidad el aprender los fundamentos del lenguaje de R en cuanto a estadística.

La razón por la que incluimos este módulo, de hecho, es porque si hay que destacar una ventaja sobre todas las demás de este lenguaje, es la capacidad de uso estadístico de todo tipo. Absolutamente de todo tipo. Está en el propio sentido de la existencia de R servir como programa de análisis estadístico, pero es que además la ingente comunidad de usuarios y desarrolladores ha permitido que las herramientas más innovadoras de la estadística estén disponibles en este programa.

En este primer documento, vamos a explorar algunas de las rutinas de código relacionadas con las distribuciones estadísticas, puesto que el uso de un determinado test requiere que algunas asunciones sobre los datos sean válidos o no.

Vamos a dejar ya cargados el paquete vcd que nos permitirá hacer algunas operaciones adicionales, tidyverse para explorar visualmente las distribuciones de dos archivos de datos que utilizaremos a lo largo de este módulo y en el ejercicio asociado.

install.packages("vcd")
library(vcd)
library(ggplot2)

fir<-read.table("fir.txt",header=TRUE)
frogs<-read.table("frogs.txt",header=TRUE)

Estas son las variables que contiene cada uno.

str(fir)
## 'data.frame':    3240 obs. of  7 variables:
##  $ VAR1       : Factor w/ 5 levels "Blackbea","Fajita",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ WAVE_NON   : Factor w/ 2 levels "n","w": 1 1 1 1 1 1 1 1 1 1 ...
##  $ TREE_NO    : Factor w/ 348 levels "1066","1068",..: 1 1 1 1 1 1 1 1 2 2 ...
##  $ año        : int  1991 1992 1993 1994 1995 1996 1997 1998 1991 1992 ...
##  $ conos      : num  0 0 0 19 0 0 0 0 0 0 ...
##  $ DBH        : num  9.4 9.4 9.4 9.4 9.4 9.4 9.4 9.4 10.6 10.6 ...
##  $ crecimiento: num  36 24 12 15 22 18 24 17 47 33 ...
str(frogs)
## 'data.frame':    48 obs. of  5 variables:
##  $ density : int  10 10 10 10 10 10 10 10 10 10 ...
##  $ pred    : Factor w/ 2 levels "no","pred": 1 1 1 1 1 1 1 1 2 2 ...
##  $ size    : Factor w/ 2 levels "big","small": 1 1 1 1 2 2 2 2 1 1 ...
##  $ surv    : int  9 10 7 10 9 9 10 9 4 9 ...
##  $ propsurv: num  0.9 1 0.7 1 0.9 0.9 1 0.9 0.4 0.9 ...

Comenzaremos con la visualización de algunas variables o columnas de nuestro dataset de abetos:

ggplot(fir,aes(x=DBH))+geom_histogram()

ggplot(fir,aes(x=conos))+geom_histogram()

ggplot(fir,aes(x=crecimiento))+geom_histogram()

ggplot(fir,aes(x=DBH))+geom_density()

ggplot(fir,aes(x=conos))+geom_density()

ggplot(fir,aes(x=crecimiento))+geom_density()

Algunas de las variables se aproximan a la forma de una campana de Gauss que nos indica la distribución normal de los datos. La asunción de normalidad es la base para que podamos usar la gran mayoría de los procedicimientos estadísticos que nos permiten validar hipótesis. Si esto asunción no se cumple, tendremos que recurrir a la estadística no paramétrica o a modelos estadísticos que no asumen esta normalidad.

El test de shapiro permite precisamente testar la normalidad de los datos de una variable. Se encuentra ya disponible por defecto en el paquete precargado stats, mediante la función shapiro.test().

shapiro.test(fir$DBH)
## 
##  Shapiro-Wilk normality test
## 
## data:  fir$DBH
## W = 0.94797, p-value < 2.2e-16
shapiro.test(fir$conos)
## 
##  Shapiro-Wilk normality test
## 
## data:  fir$conos
## W = 0.48684, p-value < 2.2e-16
shapiro.test(fir$crecimiento)
## 
##  Shapiro-Wilk normality test
## 
## data:  fir$crecimiento
## W = 0.92889, p-value < 2.2e-16

El resultado es un valor estadístico, W, y el p-valor, el cálculo estadístico que nos indica la probabilidad de que los resultados obtenidos sean un producto del azar. En este caso concreto, la probabilidad de que los datos sigan una distribución normal. En los tres casos de arriba, el valor p-value es tan bajo que rechazamos que los datos sean normales. Normalmente pondremos un umbral de p<0.05 o 0.01 por debajo del cual rechazaremos esta o hipótesis parecidas.

Cargamos también los datos de las especies de Iris para realizar también el test sobre las cuatro variables.

data("iris")
shapiro.test(iris$Sepal.Width)
## 
##  Shapiro-Wilk normality test
## 
## data:  iris$Sepal.Width
## W = 0.98492, p-value = 0.1012
shapiro.test(iris$Petal.Width)
## 
##  Shapiro-Wilk normality test
## 
## data:  iris$Petal.Width
## W = 0.90183, p-value = 1.68e-08
shapiro.test(iris$Petal.Length)
## 
##  Shapiro-Wilk normality test
## 
## data:  iris$Petal.Length
## W = 0.87627, p-value = 7.412e-10
shapiro.test(iris$Sepal.Length)
## 
##  Shapiro-Wilk normality test
## 
## data:  iris$Sepal.Length
## W = 0.97609, p-value = 0.01018

En el caso de los datos de predación sobre ranas existe la variable “surv”. La visualización de los datos indica que no siguen una distribución como una campana de Gauss. Indicaría por tanto que el test de Shapiro producirá un resultado estadísticamente significativo, es decir, se desvía de la normalidad.

ggplot(frogs, aes(x=surv))+geom_density()

shapiro.test(frogs$surv)
## 
##  Shapiro-Wilk normality test
## 
## data:  frogs$surv
## W = 0.87966, p-value = 0.0001476

En cualquier caso los datos no parecen seguir un patrón aleatorio. Los valores de supervivencia más altos tienden a ser menos frecuentes. Es posible que los datos sigan otra distribución, como la binomial negativa (para datos donde los valores menores son más frecuentes) o la distribución de Poisson, utilizada para conteos (ej. número de visitas a una flor por un polinizador).

El paquete que hemos cargado, vcd permite cotejar si una variable determinada sigue una distribución concreta, a través de la función goodfit()que necesita especificar como parámetros la variable a testar y el tipo de distribución que queremos comprobar a la que se ajusta o no. El resultado de la función en sí mismo no es excesivamente interesante: solo nos indica el desajuste entre el dato real y el dato esperable si siguiera esa distribución. Por ello, lo que salga de ella se la vamos a pasar directamente con el operador %>% de magrittr a la función summary() que identificará la clase de datos y nos entregará un resultado: un valor del test de la Chi-cuadrado y un p-valor, que nos indica que los resultados tampoco siguen una distribución binomial negativa.

goodfit(frogs$surv,type="nbinomial")%>%
  summary()
## 
##   Goodness-of-fit test for nbinomial distribution
## 
##                      X^2 df    P(> X^2)
## Likelihood Ratio 69.5072 19 1.11054e-07

Recurrimos ahora a otro de los datasets que estamos utilizando, en este caso a biocrust, que tiene un número muy grande de datos recogidos. Vamos a comprobar si los datos de la variable de temperatura media (columna Tmean) son normales, otra vez con shapiro.test(). Pero esta función tiene una limitación: solo admite menos de 5000 datos, menos de los que tiene el dataset entero. Por lo tanto vamos a limitarnos a las 4999 primeras. También vamos a evaluar visualmente la distribución de los datos:

biocrust<-read.table("biocrust.txt",header=TRUE)
shapiro.test(biocrust$Tmean[1:4999])
## 
##  Shapiro-Wilk normality test
## 
## data:  biocrust$Tmean[1:4999]
## W = 0.94836, p-value < 2.2e-16
ggplot(biocrust,aes(x=Tmean))+geom_density()

A pesar de que el test nos indica que los datos no cumplen la normalidad, visualmente sugieren otra cosa. Quizás, ¿el test es demasiado exigente? Para comprobarlo vamos a generar datos aleatorios normalmente distribuidos, en el mismo número que los datos de costras biológicas, con su misma media y su misma desviación típica. Vamos a visualizar los histogramas superpuestos de los datos reales (en color rojo) y los datos simulados:

rnorm(n=nrow(biocrust),mean=mean(biocrust$Tmean),sd=sd(biocrust$Tmean))%>% #generamos los datos con rnorm.
  as.data.frame()%>% #pasamos los datos a un dataframe.
  ggplot(aes(x=.))+geom_histogram()+ #le pasamos los datos al parémetro x de ggplot() y le añaidmos una capa de histogama
  geom_histogram(data=biocrust,aes(x=Tmean,fill="red",alpha=0.6))+ #añadimos un histograma con los datos de biocrust
  theme_bw()+ #cambiamos el tema
  theme(legend.position="NONE") #quitamos la leyenda con esta línea
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Los datos se superponen en buena medida, por lo que visualmente verificamos que los datos, al menos son casi normales, si no del todo. Vamos a ver qué pasa si le pasamos a un test de Shapiro a estos datos generados aleatoriamente:

shapiro.test(rnorm(n=5000,mean=mean(biocrust$Tmean),sd=sd(biocrust$Tmean)))
## 
##  Shapiro-Wilk normality test
## 
## data:  rnorm(n = 5000, mean = mean(biocrust$Tmean), sd = sd(biocrust$Tmean))
## W = 0.99929, p-value = 0.04266

Enocntramos un p-valor del resultado muy bajo. Según donde pongamos el umbral, podríamos decidir que ni siquiera 5000 datos son normales. La conclusión es que sí, se trata de un test excesivamente exigente.