Intervalos de confianza paramétricos con R

Dada una muestra aleatoria simple \(\xi_1, \xi_2,\ldots, \xi_n\) proveniente de una población normal, el estadístico \[\frac{\overline{\xi} - \mu}{S/\sqrt{n}} \sim t_{n-1}\] compuesto por la media muestral \( \overline\xi\), la media poblacional \( \mu\), la (cuasi) desviación típica \(S\) y el tamaño muestral \(n\) sigue una distribución \(t\) de Student con \(n-1\) grados de libertad.

Los intervalos de confianza de la media poblacional se calculan en R mediante la función t.test(), tanto el bilateral como los unilaterales.

## Data and parameters
n <- 400 # Sample size
set.seed(123) # for reproducibility
x <- rchisq(n, df = 7) # Sampling data
alpha <- 0.05 # Significance level

## Bilateral confidence interval
t.test(x, conf.level = 1 - alpha)$conf.int
## One sided confidence intervals
t.test(x, alternative = "less", conf.level = 1 - alpha)$conf.int
t.test(x, alternative = "greater", conf.level = 1 - alpha)$conf.int

[1] 6.338586 7.013779
attr(,"conf.level")
[1] 0.95
[1]     -Inf 6.959301
attr(,"conf.level")
[1] 0.95
[1] 6.393065      Inf
attr(,"conf.level")
[1] 0.95

El siguiente estadístico \[\frac{(n-1) S^2}{\sigma^2}\sim \chi^2_{n-1}\] sigue una distribución \( \chi^{2}\) con \(n-1\) grados de libertad cuando la muestra aleatoria simple proviene de una distribución normal.

Las funciones básicas de R no incluyen función alguna que calcule el intervalo de confianza para la varianza. El desarrollo de la función pivotal nos conduce a la siguiente expresión:

## Confidence interval for the variance
df <- length(x) - 1
lower <- var(x) * df / qchisq(alpha / 2, df, lower.tail = FALSE)
upper <- var(x) * df / qchisq(1 - alpha / 2, df, lower.tail = FALSE)
c(lower = lower, upper = upper)

   lower    upper 
10.31584 13.62015

Sea \( \xi_{1}, \xi_{2}, \ldots, \xi_{n}\) una muestra aleatoria simple de una distribución de Bernouilli de parámetro \(p\). Entonces \[\frac{\hat{p} -p}{\sqrt{ \frac{\hat{p}(1-\hat{p})}{n}}} \sim Normal(0,1)\] cuando \(n\) alcanza un tamaño suficientemente grande. Para valores pequeños de \(n\) se opta por el intervalo de confianza exacto.

## Data and parameters
n <- 400 # Sample size
set.seed(123) # for reproducibility
p <- 0.7 # Population proportion
x <- sample(c(0, 1), size = n, replace = TRUE, prob = c(1 - p, p)) # Sampling data
alpha <- 0.05 # Significance level

## Bilateral confidence interval
prop.test(sum(x), n = length(x), conf.level = 1 - alpha)$conf.int

## One sided confidence intervals
prop.test(sum(x), n = length(x), alternative = "less", conf.level = 1 - alpha)$conf.int
prop.test(sum(x), n = length(x), alternative = "greater", conf.level = 1 - alpha)$conf.int

[1] 0.6649956 0.7558484
attr(,"conf.level")
[1] 0.95
[1] 0.000000 0.749391
attr(,"conf.level")
[1] 0.95
[1] 0.6726582 1.0000000
attr(,"conf.level")
[1] 0.95

Para calcular el error máximo y el tamaño muestral mínimo en la estimación de proporciones en la situación más desfavorable (\(p=0.5\)), fijamos el nivel de significación y desarrollamos la función pivotal de la proporción muestral.

## Maximum error
n <- 600 # Sample size
alpha <- 0.03 # Significance level
(error <- qnorm(1 - alpha / 2) * sqrt(0.5 * 0.5 / n))

## Required sample size
error <- 0.02
alpha <- 0.03
(n <- ceiling(qnorm(1 - alpha / 2)^2 * 0.5 * 0.5 / error^2))

[1] 0.04429678
[1] 2944

Cuando observamos dos muestras aleatorias simples independientes que proceden de distribuciones normales, conviene comprobar si poseen la misma varianza (homogéneas) o no. Por defecto supondremos heterogeneidad de ambas muestras. El intervalo de confianza para la diferencia de medias se calcula como sigue:

## Data
set.seed(654)
x <- rnorm(300, 7, 3)
y <- rnorm(200, 8, 1)
alpha <- 0.025

## Confidence interval
t.test(x, y, conf.level = 1 - alpha)$conf.int
t.test(x, y, alternative = "less", conf.level = 1 - alpha)$conf.int
t.test(x, y, alternative = "greater", conf.level = 1 - alpha)$conf.int

[1] -1.7124080 -0.8731655
attr(,"conf.level")
[1] 0.975
[1]       -Inf -0.9261273
attr(,"conf.level")
[1] 0.975
[1] -1.659446       Inf
attr(,"conf.level")
[1] 0.975

Con dos muestras aleatorias simples independientes que proceden de distribuciones normales el intervalo de confianza para el cociente de varianzas se calcula mediante la distribución F de Snedecor.

## Data
set.seed(654)
x <- rnorm(300, 7, 3)
y <- rnorm(200, 8, 1)
alpha <- 0.025

var.test(x, y, conf.level = 1 - alpha)$conf.int

[1]  6.071191 10.861347
attr(,"conf.level")
[1] 0.975

Cuando disponemos de información sobre dos muestras aleatorias simples de Bernouilli independientes, se procede de la siguiente forma para calcular el intervalo de confianza de la diferencia de proporciones:

## Data
set.seed(234)
x <- sample(c(0, 1), size = 300, replace = TRUE, prob = c(0.2, 0.8)) # Sampling data
y <- sample(c(0, 1), size = 500, replace = TRUE, prob = c(0.6, 0.4)) # Sampling data
alpha <- 0.015

## Confidence interval
prop.test(c(sum(x), sum(y)), c(length(x), length(y)), conf.level = 1 - alpha)$conf.int
prop.test(c(sum(x), sum(y)), c(length(x), length(y)), alternative = "less", conf.level = 1 - alpha)$conf.int
prop.test(c(sum(x), sum(y)), c(length(x), length(y)), alternative = "greater", conf.level = 1 - alpha)$conf.int

[1] 0.3872262 0.5394404
attr(,"conf.level")
[1] 0.985
[1] -1.0000000  0.5315212
attr(,"conf.level")
[1] 0.985
[1] 0.3951455 1.0000000
attr(,"conf.level")
[1] 0.985

En los casos anteriores se conoce la distribución teórica del estimador. Normalmente se debe a que la población de referencia presenta un comportamiento normal, o bien a que el Teorema Central del Límite permite aproximar la distribución muestral a la distribución normal. En la mayoría de las situaciones basta con aplicar los anteriores intervalos de confianza para obtener resultados fiables.

Cuando se desconoce la distribución teórica subyacente del estimador, la técnica de los remuestreos (bootstrap) la aproxima mediante remuestreos a partir de la muestra original.

## Data
n <- 431 # Sample size
set.seed(123) # Reproducibility
x <- rchisq(n, 7) # Sampling data does not follow a Normal distribution.
alpha <- 0.03

Realizamos un número elevado de remuestreos con reemplazamiento (en este caso mil remuestreos). Para cada remuestreo calculamos el estimador del parámetro que deseamos inferir. Si por ejemplo, deseamos estimar la media poblacional, calculamos la media muestral. En este caso, calculemos el intervalo de confianza para la media poblacional.

meanxboot <- replicate(1000, {
    ## replace must be TRUE!
    y <- sample(x, length(x), replace = TRUE)
    mean(y)
})

## Confidence interval
quantile(meanxboot, probs = c(alpha/2, 1-alpha/2))

    1.5%    98.5% 
6.379058 7.094514

Aventurémonos a encontrar un intervalo de confianza para el percentil treinta y siete. Seguimos los pasos anteriores, pero en vez de calcular la media en cada remuestreo, determinamos el percentil treinta y siete.

quanxboot <- replicate(1000, {
    ## replace must be TRUE!
    y <- sample(x, length(x), replace = TRUE)
    quantile(y,probs=0.37)
})
quantile(quanxboot, probs = c(alpha/2, 1-alpha/2))

    1.5%    98.5% 
4.857012 5.630769

La ventaja del remuestreo reside principalmente en su sencillez. Se utiliza cuando se desconoce la distribución teórica subyacente. Aunque recuerde que «si no aprovecha la estructura subyacente, ¡deja dinero encima de la mesa!».

• George Casella, Roger L. Berger, Statistical Inference.
• Prasanna Sahoo, Probability and Mathematical Statistics.
• Carlos Gil Bellosta, Modelización estadística: si no aprovechas la estructura subyacente, dejas dinero encima de la mesa.
• David Molina Muñoz y Ana María Lara Porras, Intervalos de confianza con R.