R‎ > ‎

Estadística

Aquí veremos ejemplos que, de alguna u otra forma muestran el uso de R para resolver diversos problemas de estadística o simplemente mostrar un uso teórico de la estadística.

Los volcanes más grandes de la historia (y otras estadísticas)

publicado a la‎(s)‎ 15 abr. 2014 14:32 por Daniel Loera

Hace unos 199 años el volcán más grande de la historia registrada, el Monte Tambora, entró en erupción, arrojando una enorme cantidad de roca fundida y cenizas a la atmósfera y en la tierra circundante.

¿Cómo se cuantifica la intensidad de una erupción volcánica? La mayoría de la gente sabe acerca de la escala de Richter que cuantifica la energía liberada en un terremoto. El análogo de una erupción volcánica es el índice de explosividad volcánica (VEI).


Datos


Los datos sobre los volcanes y sus erupciones se obtuvieron del Instituto Global de Vulcanismo Smithsonian. Los datos provienen de dos componentes: uno que detalla los volcanes y otro las especificaciones de las erupciones. Ambos fueron descargados como archivos xls, que me convertí al CSV para ahorrar espacio y también para los tiempos de carga más rápidos.

Y el que escribe se los deja fácil, pues los guarde en el formato de datos de R y aquí están  para erupciones y para volcanes

> load("volcanoes.RData")
> load("eruptions.RData")
##Si no tienes los datos en el directorio de trabajo, puedes poner por ejemplo.
##load("D:/Inegi/R/tempo/volcanoes.RData")
##load("D:/Inegi/R/tempo/eruptions.RData")
> dim(volcanoes)
[1] 1550   12
> dim(eruptions)
[1] 10762    22

> table(eruptions$VEI, useNA = "ifany")

   0    1    2    3    4    5    6    7 <NA> 
 981 1271 3859 1103  492  176   51    7 2822 

Así, de los 10.762 erupciones de los datos, 2.822 de ellos no tienen un valor para VEI.

Un histograma muestra la distribución rápida de VEIs.

library(ggplot2)
ggplot(eruptions, aes(x = VEI))+geom_histogram()+ylab("")+theme_classic()

Vemos que la gran mayoría de las erupciones tienen VEI nivel 2, lo que significa que expulsan entre 0.001 y 0.01 km3 de ceniza y lava.


¿Con qué frecuencia ocurren las erupciones en función de VEI? Para obtener algunas estadísticas decentes vamos a agregar los datos por siglo. Lo primero que tenemos que hacer es agregar una columna de siglo a los datos de las erupciones. Luego contamos el número de erupciones en cada VEI para cada siglo abarcado por los datos.

> eruptions$Century = ceiling(eruptions$Start.Year / 100)
> library(plyr)
> eruption.count = ddply(eruptions, .(VEI, Century), summarise, count = length(VEI))
> eruption.count = subset(eruption.count, !is.na(VEI))
> head(eruption.count)
  VEI Century count
1   0     -95     1
2   0     -90     1
3   0     -86     1
4   0     -85     2
5   0     -84     1
6   0     -83     2
> tail(eruption.count)
    VEI Century count
527   7     -56     1
528   7     -43     1
529   7     -16     1
530   7      10     1
531   7      13     1
532   7      19     1

Una parcela de azulejos de VEI frente a los siglo da una buena vista de los datos.

ggplot(eruption.count, aes(x = Century, y = VEI, fill = log10(count)))+
  geom_tile()+scale_fill_gradientn(colours = rev(rainbow(7)),labels = 10**(0:3))+
  scale_x_continuous(breaks = seq(-100, 20, 20))+
  theme_classic()+theme(text = element_text(size = 16))


Ahora, parece que las erupciones son cada vez más frecuentes. Por supuesto, esto no es cierto: lo que sucede es que se tienen mejores registros. Los datos para el siglo pasado son bastante robusto.

> subset(eruption.count, Century==20)
    VEI Century count
95    0      20   442
124   1      20   812
208   2      20  1634
296   3      20   333
394   4      20    53
486   5      20     9
525   6      20     3
##por curiosidad puede explorar otro tipo de datos, por ejemplo.
> subset(volcanoes, Country=="Mexico")

Podemos ver que hay poco más de 3.000 erupciones con 3 o menos VEI, lo que corresponde a un promedio de alrededor de 30 por año. Erupciones con mayor VEI son sensiblemente más raro.

A pesar de la frecuencia de un tanto perturbadora de estos eventos, la mayoría de ellos ocurren en lugares remotos. Y eso se debe, probablemente, porque sabiamente decidimos no vivir demasiado cerca de ellos!

Cálculo de intervalos de confianza para proporciones

publicado a la‎(s)‎ 11 abr. 2014 11:21 por Daniel Loera

Mostraremos un par de funciones para el cálculo de los intervalos de confianza para las proporciones. En primer lugar veremos el Método simple Asymtotic :


simpasym <- function(n, p, z=1.96, cc=TRUE){
  out <- list()
  if(cc){
    out$lb <- p - z*sqrt((p*(1-p))/n) - 0.5/n
    out$ub <- p + z*sqrt((p*(1-p))/n) + 0.5/n
  } else {
    out$lb <- p - z*sqrt((p*(1-p))/n)
    out$ub <- p + z*sqrt((p*(1-p))/n)
  }
  out
}

#Usando así la función simpasym
simpasym(n=30, p=0.3, z=1.96, cc=TRUE)

Donde:
             n=    es el tamaño de la muestra
                p=    es la proporción
                z=     es el valor para el intervalo % (es decir 1.96 proporciona el IC del 95 %)
                cc=   corrección de continuidad que debe aplicarse

Los resultados devueltos son el límite inferior ($lb) y el límite superior ( $ub ) . 

El segundo método es el método de puntuación y se define de la siguiente manera:

scoreint <- function(n, p, z=1.96, cc=TRUE){
  out <- list()
  q <- 1-p
  zsq <- z^2
  denom <- (2*(n+zsq))
  if(cc){ 
    numl <- (2*n*p)+zsq-1-(z*sqrt(zsq-2-(1/n)+4*p*((n*q)+1)))
    numu <- (2*n*p)+zsq+1+(z*sqrt(zsq+2-(1/n)+4*p*((n*q)-1)))
    out$lb <- numl/denom
    out$ub <- numu/denom
    if(p==1) out$ub <- 1
    if(p==0) out$lb <- 0
  } else {
    out$lb <- ((2*n*p)+zsq-(z*sqrt(zsq+(4*n*p*q))))/denom
    out$ub <- ((2*n*p)+zsq+(z*sqrt(zsq+(4*n*p*q))))/denom
  }
  out
}

scoreint(n=40, p=0.2, z=1.96, cc=TRUE)

Veamos los resultados de simpasym scoreint con los mismos valores

simpasym(n=40, p=0.2, z=1.96, cc=TRUE)
$lb
[1] 0.06353872
$ub
[1] 0.3364613

scoreint(n=40, p=0.2, z=1.96, cc=TRUE)
$lb
[1] 0.09614404
$ub
[1] 0.3613774

Estas fórmulas ( y un par de otros) se discuten en Newcombe, RG ( 1998 ), quien sugiere que el método de puntuación debe ser con más frecuencia disponible en paquetes de software estadístico.



Referencia:
Newcombe, R. G. (1998) Two-sided confidence intervals for the single proportion: comparison of seven methods. Statist. Med., 17: 857-872. doi: 10.1002/(SICI)1097-0258(19980430)17:8<857::AID-SIM777>3.0.C

1-2 of 2