Cálculo de Intervalos de Confianza con R

• R dispone de numerosas funciones para construir intervalos de confianza con suma facilidad.

• Utilizaremos los datos del archivo datosP1Aula.xlsx (estudio del anidamiento de tortugas en Boa Vista, Cabo Verde) como ejemplo para la construcción de todos los intervalos tratados en este tema.

• Supondremos que el archivo está en la carpeta Descargas de nuestro directorio personal. Podemos leerlo mediante:

library(readxl)
tortugas <- read_excel("~/Descargas/datosP1Aula.xlsx")

Cálculo de Intervalos de Confianza con R

El contenido del archivo corresponde a datos de varias campañas para el estudio de la tortuga Caretta caretta realizadas entre 1999 y 2004 en la isla de Boa Vista (Cabo Verde). Durante cada campaña se observaron los nidos excavados en cada una de cuatro playas de la isla (Ervatao, Ponta Cosme, Calheta y Porto Ferreiro). Para cada nido se midieron las siguientes variables:

• Variables biométricas del ejemplar que construyó el nido: LCC (Longitud Curva del Caparazón, en cm.), ACC (Anchura Curva del Caparazón, en cm) y Peso.

• Tamaño de la puesta (número de huevos depositados en el nido).

• Distancia (en metros) medida perpendicularmente desde la linea de marea alta hasta la posición del nido.

• Profundidad del nido (cm).

• Para algunos nidos se midieron además: el número de crías emergidas, el número de crías muertas, el número de huevos rotos, y la presencia de cangrejos en el interior del nido (1:Sí / 0:No)

Intervalo para la media $\mu$ de una variable $X$ con distribución normal

mean(tortugas$peso)

## [1] 60.6892

t.test(tortugas$peso)$conf.int

## [1] 60.06344 61.31496
## attr(,"conf.level")
## [1] 0.95

Intervalo para la media $\mu$ de una variable $X$ con distribución normal

¿Es normal la variable peso?

La gráfica que se muestra a continuación representa un histograma de los pesos de las tortugas observadas en la isla, así como una versión "suavizada" del histograma (es lo que se denomina "densidad empírica"). Asimismo se ha superpuesto la función de densidad teórica de una variable normal cuya media y desviación típica son las correspondientes a los pesos observados. Si la variable peso es normal cabe esperar que las densidades empírica y teórica sean similares.

Intervalo para la media $\mu$ de una variable $X$ con distribución normal

ggplot(tortugas,aes(x=peso)) +
  geom_histogram(aes(y=..density..), bins=10, fill="lightcyan",
            col="black") + geom_density(aes(color="Empírica")) +  
  stat_function(fun = dnorm, aes(color="Teórica"),
                args = list(mean = mean(tortugas$peso,na.rm=TRUE), 
                            sd = sd(tortugas$peso,na.rm=TRUE)))+
  scale_color_manual("Densidad",values=c("blue","red"))

Es admisible admitir la normalidad de la variable Peso.

Intervalo para la media $\mu$ de una variable $X$ con distribución normal

Ejercicio

Intervalo para la varianza $\sigma^{2}$ de una variable con distribución normal

var(tortugas$peso)

## [1] 25.23639

library(TeachingDemos)
sigma.test(tortugas$peso, conf.level = 0.99)$conf.int

## [1] 20.25544 32.17940
## attr(,"conf.level")
## [1] 0.99

Intervalo para la varianza $\sigma^{2}$ de una variable con distribución normal

Ejercicio

Intervalo para el cociente de varianzas $\sigma_{1}^{2}/\sigma_{2}^{2}$ de variables normales.

tortugasCPC <- subset(tortugas,
                      playa=="Calheta"|playa=="Ponta Cosme") 
aggregate(peso~playa,data=tortugasCPC,var)

##         playa     peso
## 1     Calheta 35.51222
## 2 Ponta Cosme 24.04994

Cociente de varianzas $\sigma_{1}^{2}/\sigma_{2}^{2}$ de variables normales.

var.test(peso~playa, data=tortugasCPC)$conf.int

## [1] 0.868942 2.697037
## attr(,"conf.level")
## [1] 0.95

El intervalo contiene al 1. Por tanto, con una confianza del 95% no podemos asegurar que la varianza del peso sea distinta en ambas playas.

Intervalo para el cociente de varianzas $\sigma_{1}^{2}/\sigma_{2}^{2}$ de variables normales.

Ejercicio

Intervalo para la diferencia de medias de variables normales independientes

aggregate(peso~playa,data=tortugasCPC,mean)

##         playa     peso
## 1     Calheta 59.38857
## 2 Ponta Cosme 61.81667

Intervalo para la diferencia de medias de variables normales independientes

• Si suponemos que las varianzas del peso en ambas playas son iguales $\left(\sigma_{1}=\sigma_{2}\right)$, , el intervalo es:

$${\small\bbox[#ffff80]{\mu_1-\mu_2\in\left[{\left({\bar{X}_{1}-\bar{X}_{2}}\right)\pm t_{n_{1}+n_{2}-2,\alpha/2}s_{p}\sqrt{\frac{1}{{n_{1}}}+\frac{1}{{n_{2}}}}}\right]} \;\quad s_{p}=\sqrt{\frac{{\left({n_{1}-1}\right)s_{1}^{2}+\left({n_{2}-1}\right)s_{2}^{2}}}{{n_{1}+n_{2}-2}}}}$$

que en R se calcula mediante:

t.test(peso~playa,data=tortugasCPC,var.equal=TRUE)$conf.int

## [1] -4.4852781 -0.3709124
## attr(,"conf.level")
## [1] 0.95

Intervalo para la diferencia de medias de variables normales independientes

• Si suponemos que las varianzas del peso en ambas playas son distintas $\left(\sigma_{1}\ne\sigma_{2}\right)$ el intervalo es:

$${\small\bbox[#ffff80]{\mu_1-\mu_2\in\left[\left({\bar{X}_{1}-\bar{X}_{2}}\right)\pm t_{n,}{}_{\alpha/2}\sqrt{\frac{{s_{1}^{2}}}{{n_{1}}}+\frac{{s_{2}^{2}}}{{n_{2}}}}\right]} \quad \; n=\frac{{\left({\frac{{s_{1}^{2}}}{{n_{1}}}+\frac{{s_{2}^{2}}}{{n_{2}}}}\right)^{2}}}{{\left({\frac{{s_{1}^{2}}}{{n_{1}}}}\right)^{2}\frac{1}{{n_{1}-1}}+\left({\frac{{s_{2}^{2}}}{{n_{2}}}}\right)^{2}\frac{1}{{n_{2}-1}}}}}$$

que en R se calcula mediante:

t.test(peso~playa,data=tortugasCPC,var.equal=FALSE)$conf.int

## [1] -4.6991167 -0.1570738
## attr(,"conf.level")
## [1] 0.95

Este intervalo es algo más amplio que el anterior (se calcula con menos información al no asumir que las varianzas son iguales). Este intervalo es preferible al anterior

Diferencia de medias de variables no normales independientes (sólo si $\small n_{1}\ge30,n_{2}\ge30$).

Para calcular este intervalo también podemos usar la función t-test de R (considerando varianzas distintas), ya que la forma del intervalo es idéntica salvo que se ha sustituido la $t_{n-1,\alpha/2}$ por la $z_{\alpha/2}$ y para $n$ grande, ambos valores son prácticamente idénticos, $t_{n-1,\alpha/2}\cong z_{\alpha/2}$.

Diferencia de medias de variables no normales independientes (sólo si $\small n_{1}\ge30,n_{2}\ge30$).

Construimos una gráfica con los histogramas de la variable distancia en cada una de las dos playas para valorar la normalidad de dicha variable en cada playa:

ggplot(subset(tortugas,playa=="Ervatao"|playa=="Ponta Cosme"), 
       aes(x=distancia)) +
  geom_histogram(aes(y=..density..), bins=10, fill="lightcyan",
            col="black") + 
  geom_density(aes(color="Empírica")) +  
  stat_function(fun = dnorm, aes(color="Teórica"),
          args = list(mean = mean(tortugas$distancia,na.rm=TRUE), 
                      sd = sd(tortugas$distancia,na.rm=TRUE)))+
  scale_color_manual("Densidad",values=c("blue","red")) + 
  facet_grid(playa~.)

Diferencia de medias de variables no normales independientes (sólo si $\small n_{1}\ge30,n_{2}\ge30$).

La variable distancia a la orilla parece presentar una distribución no normal; no obstante, dado que en cada playa el número de observaciones es mayor que 30, podemos aplicar t-test.

Diferencia de medias de variables no normales independientes (sólo si $\small n_{1}\ge30,n_{2}\ge30$).

library(dplyr); library(knitr); library(kableExtra)
tortugas%>% 
  filter(playa=="Ervatao"|playa=="Ponta Cosme") %>% 
  group_by(playa) %>% 
  summarize(media=mean(distancia), sd=sd(distancia), n=length(distancia)) %>% 
  kable() %>%  kable_styling(full_width = FALSE, font_size=18)

t.test(distancia~playa,
       data=subset(tortugas,playa=="Ervatao"|playa=="Ponta Cosme"),
       conf.level=0.90)$conf.int

## [1] -4.083745 -0.543791
## attr(,"conf.level")
## [1] 0.9

Diferencia de medias de variables independientes

Diferencia de medias de variables normales emparejadas

Como ambas medidas se toman sobre la misma tortuga, pueden considerarse variables emparejadas. Sus valores medios en nuestra muestra de tortugas son:

mean(tortugas$LCC)

## [1] 82.2788

mean(tortugas$ACC)

## [1] 77.2448

Diferencia de medias de variables normales emparejadas

Ambas variables presentan además una correlación relativamente alta:

cor(tortugas$LCC,tortugas$ACC)

## [1] 0.8316105

ggplot(tortugas,aes(x=LCC,y=ACC)) + geom_point()

Diferencia de medias de variables normales emparejadas

$${\small\bbox[#ffff80]{\mu_{1}-\mu_{2}\in\left[\left(\overline{X}_{1}-\overline{X}_{2}\right)\pm t_{n-1,\alpha/2}\frac{S_{D}}{\sqrt{n}}\right]}\quad S_{D}=\sqrt{S_{1}^{2}+S_{2}^{2}-2S_{12}}=\sqrt{S_{1}^{2}+S_{2}^{2}-2rS_{1}S_{2}}}$$

El intervalo de confianza para la diferencia entre los valores medios de ACC y LCC en R se obtiene mediante:

t.test(tortugas$LCC, tortugas$ACC, paired=TRUE)$conf.int

## [1] 4.691706 5.376294
## attr(,"conf.level")
## [1] 0.95

Diferencia de medias de variables no normales emparejadas (sólo si $\small n\ge 60$)

Cuando $n$ es grande, los percentiles de la $t_{n-1}$ son muy parecidos a los de la $N(0,1)$:

$$z_{\alpha/2}\cong t_{n-1,\alpha/2}$$ con lo que este intervalo es prácticamente idéntico al ya visto para variables normales: $$\small \mu_{1}-\mu_{2}\in\left[\left(\overline{X}_{1}-\overline{X}_{2}\right)\pm t_{n-1,\alpha/2}\frac{S_{D}}{\sqrt{n}}\right]$$ que es el que usaremos con R en caso de variables no normales y $n$ grande.

Diferencia de medias de variables no normales emparejadas

En realidad,la variable LCC presenta indicios de falta de normalidad:

ggplot(tortugas,aes(x=LCC)) +
  geom_histogram(aes(y=..density..), bins=10) + geom_density() +  
  stat_function(fun = dnorm, color="red",
                args = list(mean = mean(tortugas$LCC,na.rm=TRUE), 
                            sd = sd(tortugas$LCC,na.rm=TRUE)))

No obstante, dado que el tamaño de muestra es grande, el intervalo calculado antes mediante t.test para la diferencia entre $\mu_{LCC}-\mu_{ACC}$ es aproximadamente correcto.

Diferencia de medias de variables no normales emparejadas (sólo si $\small n\ge 60$)

Se consideran las variables $X_1$=Distancia a la linea de marea (en cm) y $X_2=500\cdot P$ siendo $P$ la profundidad del nido (también en cm).

1. Calculamos las variables y sus medias:

X1=tortugas$distancia*1000
X2=tortugas$profNido*500
c(mean(X1),mean(X2))

## [1] 23713.6 24059.2

Diferencia de medias de variables no normales emparejadas (sólo si $\small n\ge 60$)

2. La correlación entre estas variables es:

cor(X1,X2)

## [1] 0.3693282

3. Cálculo de la diferencia entre los valores medios, así como de su intervalo de confianza:

mean(X1)-mean(X2)

## [1] -345.6

t.test(X1,X2,paired=TRUE)$conf.int

## [1] -1228.7199   537.5199
## attr(,"conf.level")
## [1] 0.95

Diferencia de medias de variables emparejadas

Intervalo de confianza para la proporción $\small\pi$ de éxitos en una $\small B(n,\pi)$

Sean $n$ el tamaño de la muestra y $N_{E}$ el número de éxitos observados.

• Método de Agresti-Coull (sólo si $n>40$, $N_E\ge5$, $n-N_E\ge5$)

$${\small\bbox[#ffff80]{\pi\in\left[\tilde{\pi}\pm z_{\alpha/2}\sqrt{\frac{\tilde{\pi}\left(1-\tilde{\pi}\right)}{\tilde{n}}}\right]},\quad \tilde{N}_{E}=N_{E}+z_{\alpha/2}^{2}/2\;\;\tilde{n}=n+\left(z_{\alpha/2}\right)^{2}\;\;\tilde{\pi}=\tilde{N}_{E}/\tilde{n}}$$

1. Comprobamos que se cumplen las condiciones de aplicación

table(tortugas$cangrejos)

## 
##   0   1 
##  72 106

En total hay $n=178$ nidos, $N_E$=106, $N-N_E=72$. Por tanto podemos aplicar el método de Agresti-Coull.

Intervalo de confianza para la proporción $\small\pi$ de éxitos en una $\small B(n,\pi)$

La proporción estimada de nidos con cangrejos es:

prop.table(table(tortugas$cangrejos))

## 
##         0         1 
## 0.4044944 0.5955056

$\hat{\pi}=0.5955$

El intervalo de Agresti-Coull es:

library(binom)
binom.confint(106,178, conf.level=0.95,method="agresti-coull")

##          method   x   n      mean     lower     upper
## 1 agresti-coull 106 178 0.5955056 0.5220969 0.6648792

Intervalo de confianza para la proporción $\small\pi$ de éxitos en una $\small B(n,\pi)$

• Método de Clopper y Pearson (cuando no se puede aplicar el anterior)

$${\small\bbox[#ffff80]{\pi\in\left[\frac{N_{E}}{{(n-N_{E}+1)F_{1}+N_{E}}},\frac{{(N_{E}+1)F_{2}}}{{(n-N_{E})+(N_{E}+1)F_{2}}}\right]},\quad\begin{array}{c} F_{1}=F_{2(n-N_{E}+1),2N_{E},\alpha/2}\\ F_{2}=F_{2(N_{E}+1),2(n-N_{E}),\alpha/2} \end{array}}$$

Calculamos primero el número de nidos con y sin cangrejos por playa:

with(tortugas, addmargins(table(playa, cangrejos)))

##                 cangrejos
## playa              0   1 Sum
##   Calheta          7  18  25
##   Ervatao         34  30  64
##   Ponta Cosme     22  47  69
##   Porto Ferreiro   9  11  20
##   Sum             72 106 178

Intervalo de confianza para la proporción $\small\pi$ de éxitos en una $\small B(n,\pi)$

El método de Clopper y Pearson se activa mediante method=exact:

binom.confint(18,25, conf.level=0.95,method="exact")

##   method  x  n mean     lower     upper
## 1  exact 18 25 0.72 0.5061232 0.8792833

binom.confint(30,64, conf.level=0.95,method="agresti-coull")

##          method  x  n    mean     lower     upper
## 1 agresti-coull 30 64 0.46875 0.3517474 0.5892916

binom.confint(47,69, conf.level=0.95,method="agresti-coull")

##          method  x  n      mean     lower    upper
## 1 agresti-coull 47 69 0.6811594 0.5637571 0.779454

binom.confint(11,20, conf.level=0.95,method="exact")

##   method  x  n mean     lower     upper
## 1  exact 11 20 0.55 0.3152781 0.7694221

Intervalo de confianza para la proporción $\small\pi$ de éxitos en una $\small B(n,\pi)$

Comparación de proporciones de dos variables binomiales $X_1$ y $X_2$

Se han realizado $n_1$ observaciones de $X_1$ con $N_{E_1}$ éxitos y $n_2$ observaciones de $X_2$ con $N_{E_2}$ éxitos. Sean $\hat{\pi_1}=\frac{N_{E_1}}{n_1}$ y $\;\hat{\pi_2}=\frac{N_{E_2}}{n_2}$ las proporciones observadas de éxitos. Si $n_{1}\ge30$, $n_{2}\ge30$:

• Intervalo de confianza para la diferencia de proporciones

  $$\left[{\left(\pi_{1}-\pi_{2}\right)\pm\left(z_{\alpha/2}\sqrt{\frac{{\hat{\pi}_{1}\left({1-\hat{\pi}_{1}}\right)}}{n_{1}}+\frac{{\hat{\pi}_{2}\left({1-\hat{\pi}_{2}}\right)}}{n_{2}}}+\frac{1}{2}\left(\frac{1}{n_{1}}+\frac{1}{n_{2}}\right)\right)}\right]$$

• Intervalo de confianza para el cociente de proporciones

  $$\left(\frac{\pi_{1}}{\pi_{2}}\right)\in\left[\left(\frac{\hat{\pi}_{1}}{\hat{\pi}_{2}}\right)\cdot\exp\left(\pm z_{\alpha/2}\sqrt{\frac{\left(1-\hat{\pi}_{1}\right)}{n_{1}\hat{\pi}_{1}}+\frac{\left(1-\hat{\pi}_{2}\right)}{n_{2}\hat{\pi}_{2}}}\right)\right]$$

Comparación de proporciones de dos variables binomiales $X_1$ y $X_2$

tb <- with(subset(tortugas,playa=="Ervatao"|playa=="Ponta Cosme"),
     table(playa,cangrejos))
tb

##              cangrejos
## playa          0  1
##   Ervatao     34 30
##   Ponta Cosme 22 47

ptb <- prop.table(tb,1)
ptb

##              cangrejos
## playa                 0         1
##   Ervatao     0.5312500 0.4687500
##   Ponta Cosme 0.3188406 0.6811594

Comparación de proporciones de dos variables binomiales $X_1$ y $X_2$

2. Su diferencia es:

diff(ptb)

##              cangrejos
## playa                  0          1
##   Ponta Cosme -0.2124094  0.2124094

3. Y el intervalo de confianza para la diferencia:

prop.test(tb)$conf.int

## [1] 0.03291732 0.39190152
## attr(,"conf.level")
## [1] 0.95

Comparación de proporciones de dos variables binomiales $X_1$ y $X_2$

1. Calculamos totales por playa en la tabla anterior:

addmargins(tb,2)

##              cangrejos
## playa          0  1 Sum
##   Ervatao     34 30  64
##   Ponta Cosme 22 47  69

2. El intervalo de confianza para cociente se obtiene mediante:

library(asbio)
ci.prat(30,64,47,69,conf=0.90,method="katz")

## 
## 90% Confidence interval for ratio of binomial proportions (method=Katz-log) 
##   Estimate        5%       95%
##  0.6881649 0.5319818 0.8902014

Comparación de proporciones de dos variables binomiales $X_1$ y $X_2$

Intervalo de confianza para el parámetro $\lambda$ de una distribución de Poisson.

siendo: $n_{1}=2T,\,\,\,\,\,n_{2}=2(T+1),\,\,\,\,\,T=\sum_{i=1}^{n}X_{i}$

Intervalo de confianza para el parámetro $\lambda$ de una distribución de Poisson.

R tampoco cuenta con ninguna función específica para el cálculo de este intervalo, pero es fácil de implementar:

intDistPoisson <- function(muestra,conf.level=0.95){
  alfa=1-conf.level
  T=sum(muestra,na.rm=TRUE)
  n=length(which(!is.na(muestra)))
  n1=2*T
  n2=2*(T+1)
  intervalo=c(qchisq(alfa/2,n1),qchisq(1-alfa/2,n2))/(2*n)
  attributes(intervalo) <- list(conf.level=conf.level)
  intervalo
}

Intervalo de confianza para el parámetro $\lambda$ de una distribución de Poisson.

Como en una distribución de Poisson de parámetro $\lambda$ se tiene que $E[X]=\lambda$, el estimador de $\lambda$ por el método de los momentos es simplemente $\hat{\lambda}=\bar{x}$.

Por tanto, en nuestro caso, el valor estimado del parámetro de la distribución del número de huevos por nido es simplemente:

mean(tortugas$Huevos)

## [1] 83.228

Calculamos su intervalo de confianza mediante la función que hemos construido antes:

intDistPoisson(tortugas$Huevos)

## [1] 82.10092 84.36668
## attr(,"conf.level")
## [1] 0.95

Intervalo de confianza para la esperanza $\mu$ de una distribución exponencial

En R no existe ninguna función para calcular este intervalo. No obstante, es muy fácil de implementar:

intDistExpo <- function(muestra,conf.level=0.95){
  media=mean(muestra,na.rm=TRUE)
  n=length(which(!is.na(muestra)))
  alfa=1-conf.level
  intervalo=2*n*media/qchisq(c(1-alfa/2,alfa/2),2*n)
  attributes(intervalo) <- list(conf.level=conf.level)
  intervalo
}

Intervalo de confianza para la esperanza $\mu$ de una distribución exponencial

En nuestra muestra de tortugas no hay ninguna variable cuya distribución sea exponencial. Por ello, simulamos a continuación 100 valores de una variable exponencial de parámetro $\mu=10$ y mostramos los 50 primeros valores resultantes:

x=rexp(100,1/10)
x[1:50]

##  [1]  4.3788474 29.2095174  1.2351388 25.1886679  2.8538529 16.6726073
##  [7]  4.1436835  9.3326899  4.3362465  1.8705817  4.3309157  1.3438427
## [13]  5.9837254 40.4264501  3.5569852  2.2283670 13.8276553  8.6614246
## [19]  7.2104161  1.2678965 21.7669334  0.3983911 12.2596748  2.1979290
## [25] 37.6819227  5.5277873  0.4728925  1.2937034  6.0308080 28.4465475
## [31]  2.9329699  1.8517686  0.6204385  1.0672398  8.4053534  2.4600792
## [37] 13.8599722 17.7843863 11.5689159  7.7059711 19.3221310  2.0510295
## [43]  9.4370670  0.7178946 17.8472992  6.1104391  0.3973279 13.0910297
## [49]  0.3169744  2.1066865

Intervalo de confianza para la esperanza $\mu$ de una distribución exponencial

Dibujamos el histograma de la variable anterior; efectivamente su aspecto es de una exponencial negativa:

ggplot(data.frame(x=x), aes(x)) + geom_histogram(bins=10)

Intervalo de confianza para la esperanza $\mu$ de una distribución exponencial

Calculamos la media observada de los valores de esta variable (recuérdese que en el caso de la exponencial $E[X]=\mu$, y por tanto podemos estimar $\hat{\mu}=\bar{x}$):

mean(x)

## [1] 9.367754

Por último utilizamos nuestra función para obtener el intervalo de confianza:

intDistExpo(x)

## [1]  7.772202 11.513390
## attr(,"conf.level")
## [1] 0.95