01 / Fundamentos

De la pregunta científica a la decisión estadística

La población es el conjunto de unidades o el proceso al que se refiere la investigación; la muestra es la colección de observaciones disponibles. Un parámetro, como la media poblacional μ o la varianza poblacional σ², caracteriza la población. Un estimador, como la media muestral X̄ o la varianza muestral S², es una función de la muestra. Su valor calculado con los datos concretos es una estimación.

La distinción es esencial: los datos son observados; el parámetro suele ser desconocido; el estimador varía de una muestra a otra. La inferencia estudia esa variación bajo un modelo y un diseño de obtención de datos. Una muestra grande no corrige por sí sola un sesgo de selección ni vuelve independientes observaciones que no lo son.

La desviación estándar resume la dispersión de las observaciones. El error estándar describe la dispersión de un estimador entre muestras. Para observaciones independientes e idénticamente distribuidas con varianza finita, el error estándar de la media es σ/√n; al desconocerse σ, se estima mediante S/√n. Esta expresión necesita ajustes cuando hay dependencia, conglomerados, ponderaciones o un diseño muestral complejo.

Hipótesis, estadístico y valor p

Una hipótesis nula H₀ delimita el modelo de referencia. Para una media puede formularse H₀: μ = μ₀. La alternativa bilateral es H₁: μ ≠ μ₀; una alternativa unilateral especifica μ > μ₀ o μ < μ₀. La dirección debe responder a la pregunta científica y establecerse antes de utilizar el resultado observado para elegirla.

El estadístico resume la discrepancia entre los datos y H₀ en una escala cuya distribución puede estudiarse. El valor p es, bajo H₀ y los demás supuestos del modelo, la probabilidad de obtener un estadístico al menos tan extremo como el observado, según la alternativa adoptada. No es la probabilidad de que H₀ sea verdadera ni la probabilidad de que el resultado se deba al azar. La significación tampoco mide el tamaño o la importancia del efecto.^[1]

El nivel α fija un límite de error tipo I de la regla; en procedimientos conservadores, la probabilidad efectiva puede ser menor. β depende del efecto verdadero, del tamaño muestral, de la variabilidad y de la prueba. No rechazar H₀ puede expresar compatibilidad con ella, falta de precisión o poca capacidad para detectar un efecto relevante. No demuestra igualdad.

El resultado estadístico adquiere significado al volver al fenómeno, a la unidad de medida y a las condiciones en que se produjeron los datos.

02 / Preparación

Datos, importación y estructura

El ejemplo original considera una muestra aleatoria de 21 tiendas de diferentes partes del país. Registra un indicador de limpieza, un índice de producto reparado y comercializado, y las unidades vendidas en el último mes. En esta edición se conserva ese escenario como ejercicio; los valores siguientes son ilustrativos. No se infiere de ellos la situación de una población real.

En un proyecto de RStudio conviene trabajar con rutas relativas al proyecto. getwd() informa el directorio actual; setwd("ruta/real") permite cambiarlo. La cadena "(...)" del texto original es un marcador, no una ruta ejecutable. La siguiente alternativa de importación presupone que el archivo no contiene encabezado y tiene tres columnas separadas por espacios o tabulaciones:

getwd()
# setwd("ruta/real/al/proyecto")  # solo si necesita cambiarla
# Si dispone del archivo original, sin encabezado:
# estudios <- read.table("estudios.txt", header = FALSE,
#                        col.names = c("menor16", "ipc", "ventas"))

read.table() lee el archivo; col.names asigna nombres a las columnas del marco de datos. No crea el archivo físico. Si este ya contiene nombres, corresponde usar header = TRUE; si tiene otro separador o marca decimal, deben ajustarse sep y dec. Leerlo dos veces, primero sin asignarlo y luego asignándolo, no es necesario.

Para ejecutar el artículo sin ese archivo, se utiliza el conjunto completo siguiente:

estudios <- data.frame(
    menor16 = c(12,15,11,18,14,16,13,19,17,12,20,15,18,14,16,19,13,17,11,20,18),
    ipc = c(11,12,14,16,13,15,10,16,14,12,18,20,17,22,19,24,18,21,23,17,20),
    ventas = c(145,158,172,183,167,192,151,176,164,180,188,205,194,217,201,229,182,212,221,196,208)
  )

Los comandos names(), nrow(), ncol(), dim() y str() permiten verificar nombres, número de filas, número de columnas, dimensiones y estructura. Aquí deben obtenerse 21 filas y tres columnas numéricas.

stopifnot(all(c("menor16", "ipc", "ventas") %in% names(estudios)),
          all(vapply(estudios, is.numeric, logical(1))))
print(names(estudios)); print(nrow(estudios)); print(ncol(estudios))
print(dim(estudios)); str(estudios)
print(colSums(is.na(estudios)))
if (any(!is.finite(as.matrix(estudios)))) {
  stop("Hay valores ausentes o no finitos: documente su tratamiento antes de analizar.")
}
stopifnot(nrow(estudios) >= 2L)

Acceso explícito y el recurso histórico a attach()

attach(estudios) incorpora las variables a la ruta de búsqueda y permite escribir sus nombres sin el prefijo del marco de datos; detach(estudios) retira esa entrada. Es útil entender esos comandos porque aparecen en el original, pero un objeto homónimo puede ocultar otro y las modificaciones posteriores del marco no se reflejan necesariamente como se espera. En el código revisado se utiliza estudios$ventas, with(estudios, ...) o el argumento data = estudios. La asignación ventas <- estudios$ventas crea un vector explícito y debe actualizarse si luego se modifica el marco.

03 / Exploración

Medias, dispersión, cuantiles e histogramas

La descripción precede al contraste. La media y la desviación estándar informan aspectos distintos; los cuantiles y las representaciones gráficas permiten detectar asimetrías, extremos o agrupamientos que un único resumen puede ocultar.

ventas <- estudios$ventas
mean(ventas)
sd(ventas)
var(ventas)
apply(estudios, 2, mean)
apply(estudios, 2, sd)
quantile(ventas)
apply(estudios, 2, quantile)
qv <- quantile(ventas, probs = c(0.025, 0.975))
qv

apply(estudios, 2, mean) aplica la función a cada columna; MARGIN = 1 la aplicaría por filas. No es necesario ejecutar previamente mean() o quantile(). En un marco de datos mixto, apply() puede convertir todo a una matriz de caracteres; por eso aquí se verifica que las columnas sean numéricas. vapply() y lapply() permiten trabajar por columnas sin esa conversión global.^[2]

Qué expresan los cuantiles

Los cuantiles localizan valores de la variable a partir de proporciones acumuladas. quantile(ventas) devuelve por defecto mínimo, cuartiles y máximo; probs = c(0.025, 0.975) solicita los percentiles 2,5 y 97,5. R dispone de varios algoritmos; el predeterminado es type = 7, basado en interpolación.^[3]

En los registros ilustrativos, la media es 187,667, la desviación estándar 23,258 y esos dos percentiles son 148 y 225. Los valores 141 y 243 mencionados en el original no se reproducen aquí porque corresponden a otro archivo de datos. Con muestras pequeñas, los cuantiles de cola son especialmente sensibles a pocas observaciones.

tinta <- "#2c2418"; cobre <- "#8b4513"; beige <- "#d4c8b0"
par(col = tinta, col.axis = tinta, col.lab = tinta, col.main = tinta,
    bg = "#fffdf7", las = 1)
hist(ventas, breaks = 7, col = beige, border = cobre,
     xlab = "Ventas", ylab = "Frecuencia", main = "Ventas: datos ilustrativos")
abline(v = qv, col = cobre, lty = 2)

abline(v = qv) dibuja líneas verticales en los cuantiles; col modifica el color y lty el tipo de línea. Un valor numérico de breaks propone un número de clases, pero R puede elegir límites convenientes y devolver un número diferente. Para exigir límites precisos, se proporciona un vector explícito de puntos de corte. Si las clases tienen anchos diferentes, es importante distinguir frecuencias de densidades.

04 / Simulación

Poblaciones normales y funciones propias

Adoptamos la notación X ∼ N(μ, σ²), donde el segundo parámetro es la varianza. En R, sin embargo, rnorm(n, mean, sd) recibe la desviación estándar. Por tanto, rnorm(1000, 35, 2) genera 1 000 observaciones de N(35, 4). En esta llamada 35 y 2 son parámetros de la distribución generadora, no la media y desviación estándar efectivamente observadas en la muestra.^[4]

set.seed(20210906)
y <- rnorm(1000, mean = 35, sd = 2)
qy <- quantile(y, probs = c(0.025, 0.975))
hist(y, probability = TRUE, col = beige, border = cobre,
     main = "N(35, 4): densidad estimada y teorica", xlab = "Valor")
abline(v = qy, col = cobre, lty = 2)
lines(density(y), col = cobre, lwd = 2)
curve(dnorm(x, mean = 35, sd = 2), add = TRUE, lty = 3, lwd = 2)
legend("topright", c("Estimada por nucleos", "Normal teorica"),
       col = c(cobre, tinta), lty = c(1, 3), bty = "n", cex = 0.7)

La semilla de set.seed() permite repetir la secuencia pseudoaleatoria dentro de un entorno reproducible. Las versiones de R y la configuración del generador también importan; sessionInfo() documenta la sesión. El histograma con probability = TRUE expresa densidades y tiene área total unitaria, lo que permite compararlo con curvas de densidad.

density(y) estima una curva a partir de los datos mediante núcleos y un ancho de banda. No conoce la densidad normal que generó la muestra. Para representar esa densidad teórica se utiliza dnorm(). La coincidencia aproximada entre ambas es un resultado de la simulación, no una identidad.^[5]

Una función para simular, graficar y devolver resultados

La función plot.m() conserva el propósito del ejemplo original, agrega una línea para la media, comprueba sus entradas y devuelve los datos y resúmenes. invisible() evita una impresión extensa automática, pero permite guardar el resultado.

plot.m <- function(n, m, s) {
  stopifnot(length(n) == 1L, is.finite(n), n >= 2, n == as.integer(n),
            length(m) == 1L, is.finite(m), length(s) == 1L, is.finite(s), s > 0)
  y <- rnorm(n, mean = m, sd = s)
  qy <- quantile(y, c(0.025, 0.975))
  hist(y, probability = TRUE, col = beige, border = cobre,
       main = paste("Muestra normal: n =", n), xlab = "Valor")
  abline(v = qy, col = cobre, lty = 2)
  abline(v = mean(y), col = tinta, lwd = 2)
  lines(density(y), col = cobre)
  invisible(list(media = mean(y), sd = sd(y), cuantiles = qy, datos = y))
}
muestra_10 <- plot.m(10, 100, 15)
muestra_10000 <- plot.m(10000, 100, 15)
print(muestra_10[c("media", "sd", "cuantiles")])
print(muestra_10000[c("media", "sd", "cuantiles")])

En la versión original un encabezado anunciaba n = 10, mientras la llamada utilizaba 10 000. Aquí se muestran ambos tamaños, con el mismo modelo N(100, 15²), como experiencias distintas. Una muestra mayor suele producir una descripción más estable de la población; no hace que las observaciones individuales tengan menor dispersión poblacional.

Las cuatro familias de funciones de distribución merecen distinguirse: d* evalúa una densidad o masa, p* una probabilidad acumulada, q* un cuantil y r* genera realizaciones. Así, dnorm(), pnorm(), qnorm() y rnorm() son operaciones diferentes sobre una misma familia de distribuciones.

05 / Variación entre muestras

Distribución muestral y error estándar de la media

Una población, una muestra de ella y una distribución de medias muestrales son objetos distintos. Para estudiar el tercero se repite un mismo procedimiento: tomar n observaciones, calcular su media, guardar el resultado y comenzar con otra muestra independiente.

En este caso la distribución normal de la media es exacta para cualquier n. No hace falta invocar una aproximación asintótica para obtenerla. Con una población no normal, el teorema central del límite establece una aproximación bajo sus condiciones; no existe un umbral universal de tamaño muestral que garantice una buena aproximación en todos los fenómenos.

n_sim <- 10000L; mu_sim <- 100; sigma_sim <- 15; B <- 1000L
medias <- numeric(B)
for (i in seq_len(B)) {
  sam <- rnorm(n_sim, mean = mu_sim, sd = sigma_sim)
  medias[i] <- mean(sam)
}
print(c(media_teorica = mu_sim, media_simulada = mean(medias),
        ee_teorico = sigma_sim / sqrt(n_sim), ee_simulado = sd(medias)))
qm <- quantile(medias, c(0.025, 0.975))
hist(medias, probability = TRUE, col = beige, border = cobre,
     main = "Distribucion de medias", xlab = "Media muestral")
abline(v = qm, col = cobre, lty = 2); lines(density(medias), col = cobre)

Con σ = 15 y n = 10 000, el error estándar teórico es 0,15. sd(medias) lo aproxima por simulación; no tiene por qué coincidir exactamente porque se han usado 1 000 repeticiones, no infinitas. La media de esas 1 000 medias debe quedar cerca de 100. Su error Monte Carlo teórico es 0,15/√1000, aproximadamente 0,00474.

numeric(B) reserva espacio para B resultados. En el bucle for, sam contiene la muestra completa de la repetición actual; medias[i] guarda solamente su media. El bucle repite instrucciones: no es por sí mismo una función, aunque puede escribirse dentro de una. La analogía con órbitas computacionales puede ayudar a representar la reiteración, pero no implica que el procedimiento sea un modelo dinámico del fenómeno investigado.

Extensión: una población inicialmente asimétrica

Una distribución exponencial de tasa 1 tiene media 1 y varianza 1. Comparar las medias de muestras de tamaños 1, 5 y 30 permite observar cómo cambia su forma y cómo su desviación estándar se aproxima a 1/√n. Este ejercicio ilustra un proceso de aproximación distinto de la normalidad exacta del ejemplo anterior.

set.seed(20210907)
medias_exp <- lapply(c(1L, 5L, 30L), function(nn) {
  replicate(4000L, mean(rexp(nn, rate = 1)))
})
names(medias_exp) <- c("n1", "n5", "n30")
print(vapply(medias_exp, mean, numeric(1)))
print(vapply(medias_exp, sd, numeric(1)))

En replicate() se repite una expresión; lapply() organiza los resultados para varios tamaños. La simulación representa matemáticamente muestreos repetidos. No sustituye un experimento real ni acredita por sí misma que un modelo describa adecuadamente observaciones de laboratorio, tiendas o pacientes.

06 / Estimación

Intervalos de confianza para una media

Para una muestra independiente de una población normal con σ desconocida, el estadístico (X̄ − μ)/(S/√n) sigue una distribución t con n − 1 grados de libertad. Al invertir su región central se obtiene el intervalo bilateral:

La distribución t surge de combinar una media normal estandarizada y una estimación independiente de su varianza. Sus colas son más pesadas que las de la normal estándar para grados de libertad finitos; al aumentar estos, se aproxima a ella. qt(0.975, n - 1) devuelve el cuantil que deja el 97,5 % de probabilidad a su izquierda; para un intervalo bilateral del 95 % quedan 2,5 % en cada cola.^[6]

n_ventas <- length(ventas)
t_critico <- qt(0.975, df = n_ventas - 1)
ee_ventas <- sd(ventas) / sqrt(n_ventas)
mean(ventas) - t_critico * ee_ventas
mean(ventas) + t_critico * ee_ventas
mean(ventas) + c(-1, 1) * t_critico * ee_ventas

El vector c(-1, 1) construye simultáneamente los dos extremos. La siguiente función generaliza el cálculo para un nivel de confianza p; este argumento no es el valor p de una prueba. Admite la exclusión explícita de ausencias, exige al menos dos observaciones y detiene el cálculo si no hay variación.

ic <- function(x, p = 0.95, na.rm = FALSE) {
  stopifnot(is.numeric(x), length(p) == 1L, is.finite(p), p > 0, p < 1,
            is.logical(na.rm), length(na.rm) == 1L, !is.na(na.rm))
  if (na.rm) x <- x[!is.na(x)]
  if (any(!is.finite(x))) stop("Datos ausentes o no finitos.")
  n_x <- length(x)
  if (n_x < 2L) stop("Se requieren al menos dos observaciones.")
  if (sd(x) == 0) stop("No hay variacion muestral.")
  ee_x <- sd(x) / sqrt(n_x)
  critico <- qt((1 + p) / 2, df = n_x - 1)
  setNames(mean(x) + c(-1, 1) * critico * ee_x, c("inferior", "superior"))
}
print(ic(ventas, 0.95)); print(ic(ventas, 0.99))
print(ic(estudios$ipc, 0.95)); print(ic(estudios$menor16, 0.95))
stopifnot(isTRUE(all.equal(unname(ic(ventas)),
                          as.numeric(t.test(ventas)$conf.int))))

En los datos ilustrativos, el IC del 95 % de la media de ventas es [177,080; 198,254], mientras que el del 99 % es [173,226; 202,108]. Al mantener los datos y aumentar la confianza, aumenta el valor crítico y el intervalo se ensancha. Esto no significa que el estimador sea más robusto frente a valores extremos o errores de especificación.

La interpretación frecuentista del 95 % se refiere al procedimiento: si se repitiera el muestreo y la construcción del intervalo bajo las mismas condiciones, el 95 % de los intervalos contendría μ. Una vez observados los datos, el parámetro es fijo y el intervalo concreto lo contiene o no. Tampoco se afirma que el 95 % de las tiendas tenga ventas dentro del intervalo para μ.

07 / Un parámetro

Prueba t para la media de una muestra

La pregunta del ejercicio es si la media poblacional de ventas difiere de 180 unidades mensuales. Adoptamos H₀: μ = 180, H₁: μ ≠ 180 y α = 0,05. El valor de referencia debe tener una justificación externa —por ejemplo, un estándar o una hipótesis previa— y no elegirse para obtener el resultado deseado.

alpha <- 0.05; mu0 <- 180
n_ventas <- length(ventas); gl <- n_ventas - 1L
ee_ventas <- sd(ventas) / sqrt(n_ventas)
t_obs <- (mean(ventas) - mu0) / ee_ventas
p_manual <- 2 * pt(-abs(t_obs), df = gl)
prueba_una <- t.test(ventas, mu = mu0, alternative = "two.sided",
                    conf.level = 1 - alpha)
print(prueba_una)
print(c(t_manual = t_obs, p_manual = p_manual))
stopifnot(isTRUE(all.equal(p_manual, prueba_una$p.value)))
# Una alternativa unilateral se elige antes de mirar los resultados.
prueba_mayor <- t.test(ventas, mu = mu0, alternative = "greater")

La fórmula con abs() es válida para ambos signos del estadístico. pt() evalúa la función de distribución acumulada de la t, no una densidad. La duplicación se justifica por la simetría de esa distribución y por la alternativa bilateral. Utilizar pt(-abs(t_obs), ...) evita además la resta 1 - pt(...), que puede perder precisión en una cola muy pequeña.^[6]

Pruebas unilaterales y dualidad con los intervalos

alternative = "greater" contrasta una media mayor que μ₀; "less", una menor. Cambian la región de rechazo y el intervalo que R devuelve. No debe elegirse una cola después de observar el signo del efecto. Para el contraste bilateral y el intervalo bilateral correspondientes, al mismo nivel y bajo el mismo procedimiento, excluir μ₀ equivale a rechazar H₀, salvo la convención en la frontera.

La prueba no se realiza «dentro» de un intervalo como operación independiente: contraste e intervalo son dos formas relacionadas de resumir la misma información bajo el modelo. El valor p no reemplaza el intervalo ni la magnitud en unidades originales.

08 / Comparación

Construir grupos sin perder de vista la variable original

El ejercicio divide las tiendas según ipc < 17 e ipc ≥ 17. El umbral coherente con el código original es 17; la mención a 117 era una inconsistencia de su encabezado. Crear dos categorías no convierte el estudio en un experimento ni vuelve causal una asociación.

umbral <- 17
estudios$ipc1 <- ifelse(estudios$ipc < umbral, 1L, 2L)
estudios$ipc2 <- factor(estudios$ipc1, levels = c(1, 2),
                        labels = c("ipc < 17", "ipc >= 17"))

med_grupos <- tapply(estudios$ventas, estudios$ipc2, mean)
sd_grupos <- tapply(estudios$ventas, estudios$ipc2, sd)
var_grupos <- tapply(estudios$ventas, estudios$ipc2, var)
n_grupos <- table(estudios$ipc2)
stopifnot(length(n_grupos) == 2L, all(n_grupos >= 2L), all(var_grupos > 0))
print(cbind(n = as.numeric(n_grupos), media = med_grupos,
            sd = sd_grupos, varianza = var_grupos))

ifelse() expresa la clasificación; factor() asigna niveles y etiquetas. Una variable de agrupación clasifica observaciones, pero no es sinónimo de variable en general. tapply() aplica una función a los grupos definidos por niveles de factores; la referencia técnica a una ragged array no significa que un marco de datos tenga filas de longitudes diferentes.^[7]

En los registros construidos, el grupo con ipc < 17 tiene 10 tiendas y media 168,8; el grupo con ipc ≥ 17 tiene 11 tiendas y media 204,818. El orden de los niveles fija el sentido de la diferencia: primero menos segundo.

Diagramas de caja y medias con intervalos

El diagrama de caja resume mediana y cuartiles; sus bigotes dependen de la convención utilizada y no son intervalos de confianza de la media. Los puntos con barras de confianza representan otro objeto: estimaciones de medias e incertidumbre. gplots::plotmeans() ofrece este segundo recurso y se conserva como opción. Un marcador romboidal no convierte todo gráfico de intervalos en un «gráfico de diamantes» con un significado estadístico universal.

plot(ipc ~ ipc2, data = estudios, col = beige, border = cobre,
     xlab = "Grupo", ylab = "IPC", main = "Clasificacion por IPC")
abline(h = umbral, col = cobre, lty = 2)
boxplot(ventas ~ ipc2, data = estudios, col = beige, border = cobre,
        ylab = "Ventas", xlab = "Grupo", main = "Ventas por grupo")

ic_grupos <- sapply(split(estudios$ventas, estudios$ipc2), ic)
plot(seq_along(med_grupos), med_grupos, pch = 18, xaxt = "n",
     xlim = c(0.5, 2.5), ylim = range(ic_grupos), xlab = "Grupo", ylab = "Ventas",
     main = "Medias e intervalos de confianza del 95%")
axis(1, at = 1:2, labels = names(med_grupos))
arrows(1:2, ic_grupos[1, ], 1:2, ic_grupos[2, ],
       angle = 90, code = 3, length = 0.05, col = cobre)
# Recurso opcional del articulo original:
if (requireNamespace("gplots", quietly = TRUE)) {
  gplots::plotmeans(ventas ~ ipc2, data = estudios, connect = FALSE,
                    col = cobre, barcol = cobre)
}

La superposición o ausencia de superposición entre dos intervalos individuales del 95 % no sustituye un contraste de la diferencia entre medias. Para esa pregunta se necesita el error estándar de la diferencia y su intervalo propio. Una representación útil puede acompañar al contraste, pero no lo reemplaza mecánicamente.

09 / Dos medias

Student con varianza combinada y prueba de Welch

Se contrasta H₀: μ₁ − μ₂ = 0 frente a una alternativa bilateral. Ambas versiones requieren independencia entre las unidades; difieren en el tratamiento de la variabilidad de los grupos. Se debe mantener el orden de la resta al interpretar el signo.

Varianzas poblacionales iguales

Si se justifica σ₁² = σ₂², la varianza combinada pondera las varianzas muestrales por sus grados de libertad. No es una media aritmética simple ni es «más precisa» sin condición: esa justificación depende del modelo de varianza común.

dif_medias <- unname(med_grupos[1] - med_grupos[2])
sp2 <- sum((as.numeric(n_grupos) - 1) * var_grupos) / (sum(n_grupos) - 2)
ee_pool <- sqrt(sp2 * sum(1 / as.numeric(n_grupos)))
t_pool <- dif_medias / ee_pool
gl_pool <- sum(n_grupos) - 2
p_pool <- 2 * pt(-abs(t_pool), df = gl_pool)
prueba_pool <- t.test(ventas ~ ipc2, data = estudios, var.equal = TRUE)
print(prueba_pool)

Varianzas no impuestas como iguales: Welch

Welch estima por separado las varianzas y utiliza grados de libertad de Welch–Satterthwaite. Estos no tienen por qué ser enteros. Es el comportamiento predeterminado de t.test() para dos muestras independientes; var.equal = TRUE solicita expresamente la versión anterior.^[8]

a_w <- var_grupos / as.numeric(n_grupos)
ee_welch <- sqrt(sum(a_w))
gl_welch <- sum(a_w)^2 / sum(a_w^2 / (as.numeric(n_grupos) - 1))
t_welch <- dif_medias / ee_welch
p_welch <- 2 * pt(-abs(t_welch), df = gl_welch)
prueba_welch <- t.test(ventas ~ ipc2, data = estudios, var.equal = FALSE)
print(prueba_pool); print(prueba_welch)
stopifnot(isTRUE(all.equal(unname(p_pool), prueba_pool$p.value)),
          isTRUE(all.equal(unname(p_welch), prueba_welch$p.value)),
          isTRUE(all.equal(unname(t.test(ventas ~ ipc2, data = estudios)$p.value),
                           prueba_welch$p.value)))

El intervalo de la diferencia se construye como (x̄₁ − x̄₂) ± valor crítico × error estándar, usando los grados de libertad pertinentes. Si el intervalo bilateral correspondiente contiene cero, no se rechaza la igualdad al mismo nivel. Esto no equivale a una demostración de equivalencia: esa pregunta necesita márgenes sustantivos y un procedimiento diseñado para contrastarlos.

Una prueba F preliminar no debe funcionar como interruptor automático entre Student y Welch. No rechazar igualdad de varianzas no acredita que sean iguales; la estrategia completa de selección cambia las propiedades del procedimiento. Para este cuaderno, Welch es la opción principal y la versión combinada permite estudiar el caso de varianza común.

10 / Dos varianzas

La prueba F y el intervalo de una razón

Comparar dispersiones es una pregunta diferente de comparar medias. La hipótesis general es H₀: σ₁²/σ₂² = r₀, con r₀ > 0. Para igualdad se toma r₀ = 1. Esto no exige σ₁² = σ₂² = 1, sino solamente igualdad de las dos varianzas.

La distribución F es positiva y asimétrica. Para el contraste bilateral de colas iguales implementado aquí, el valor p es el doble de la menor probabilidad de cola, acotado por uno. Duplicar siempre pf(F_obs, ...) puede dar una probabilidad errónea e incluso mayor que uno.^[9]

razon0 <- 1
razon_muestral <- unname(var_grupos[1] / var_grupos[2])
f_obs <- razon_muestral / razon0
gl1 <- as.numeric(n_grupos[1]) - 1; gl2 <- as.numeric(n_grupos[2]) - 1
p_f <- min(1, 2 * min(pf(f_obs, gl1, gl2),
                       pf(f_obs, gl1, gl2, lower.tail = FALSE)))
ic_razon <- razon_muestral / qf(c(1 - alpha/2, alpha/2), gl1, gl2)
prueba_f <- var.test(ventas ~ ipc2, data = estudios, ratio = razon0)
print(prueba_f); print(ic_razon)
stopifnot(isTRUE(all.equal(p_f, prueba_f$p.value)),
          isTRUE(all.equal(ic_razon, as.numeric(prueba_f$conf.int))))

El intervalo estima la razón poblacional; no está centrado simétricamente en r̂ en la escala original. Si contiene uno, no se rechaza la igualdad mediante el contraste bilateral correspondiente. Para probar una razón diferente se compara con r₀, no necesariamente con uno.

En los datos ilustrativos, la razón muestral es aproximadamente 1,049 y p ≈ 0,9338. No se rechaza la igualdad de varianzas. No se ha demostrado su identidad ni que las medias sean iguales.

11 / Dependencia dentro de pares

Dos tratamientos de sueño en las mismas personas

El conjunto histórico sleep contiene 20 observaciones de diez pacientes, con dos tratamientos soporíferos. extra es el incremento de horas de sueño respecto al control; group identifica el tratamiento y ID, al paciente. No son veinte personas independientes.^[11]

El emparejamiento debe construirse por identificador, no suponiendo que dos vectores mantienen el mismo orden. En sleep$group == "1", el símbolo $ accede a la columna group; el «1» selecciona un nivel de tratamiento y no significa «columna número uno». Las filas y columnas seleccionadas se distinguen con la coma de [filas, columnas].

data("sleep", package = "datasets")
A <- sleep[sleep$group == "1", c("ID", "extra")]
B_sleep <- sleep[sleep$group == "2", c("ID", "extra")]
stopifnot(!anyDuplicated(A$ID), !anyDuplicated(B_sleep$ID),
          setequal(A$ID, B_sleep$ID))
pares <- merge(A, B_sleep, by = "ID", suffixes = c("_A", "_B"), sort = TRUE)
stopifnot(nrow(pares) == 10L, all(complete.cases(pares)))

Definimos Dᵢ = Aᵢ − Bᵢ. La prueba pareada de H₀: μ_D = 0 es una prueba t de una muestra aplicada a las diferencias:

La covarianza explica por qué el análisis pareado puede tener una precisión diferente del análisis que ignora las parejas. Su ventaja depende de la estructura de dependencia. Las diferencias de pacientes distintos deben ser independientes; para la t exacta, la normalidad se refiere a las diferencias, no a imponer dos distribuciones marginales normales por separado.

dif_sueno <- pares$extra_A - pares$extra_B
prueba_pareada <- t.test(pares$extra_A, pares$extra_B, paired = TRUE)
prueba_diferencias <- t.test(dif_sueno, mu = 0)
prueba_independiente <- t.test(pares$extra_A, pares$extra_B)
print(prueba_pareada); print(prueba_independiente)
print(c(dz = mean(dif_sueno) / sd(dif_sueno)))
stopifnot(isTRUE(all.equal(prueba_pareada$p.value, prueba_diferencias$p.value)))
matplot(1:2, t(as.matrix(pares[c("extra_A", "extra_B")])),
        type = "b", pch = c(1, 16), lty = 1, col = cobre,
        xaxt = "n", xlab = "Tratamiento", ylab = "Incremento de horas de sueno",
        main = "Mediciones pareadas por paciente")
axis(1, at = 1:2, labels = c("A", "B"))
# Grafico por paciente conservado del ejemplo original:
plot(seq_len(nrow(pares)), pares$extra_A, type = "b", pch = 1,
     col = cobre, ylim = range(pares[c("extra_A", "extra_B")]),
     xlab = "Paciente", ylab = "Incremento de horas", main = "A y B por paciente")
lines(seq_len(nrow(pares)), pares$extra_B, type = "b", pch = 16, lty = 2)
legend("topleft", c("A", "B"), col = c(cobre, tinta), pch = c(1, 16),
       lty = c(1, 2), bty = "n")

La llamada sin paired = TRUE, conservada como comparación didáctica, produce p ≈ 0,07939 porque trata los vectores como muestras independientes. Para este diseño, esa llamada ignora información del emparejamiento. No se elige el análisis por cuál cruza 0,05: se elige por la estructura de los datos.

La diferencia estandarizada pareada d_z = d̄/s_D es aproximadamente −1,285. Esta definición no es la misma que dividir una diferencia de grupos independientes por una desviación estándar combinada; al informar un efecto estandarizado, debe identificarse su denominador.

Se conserva también el gráfico del original que compara A y B a lo largo del índice de pacientes. Ese eje representa identificadores, no una serie temporal. Un diagrama con una línea por paciente entre tratamientos hace visible el emparejamiento con mayor claridad. Los intervalos marginales de plotmeans(extra ~ group) no son el intervalo de confianza de la diferencia pareada.

12 / Diseño de investigación

Potencia, diferencia detectable y tamaño muestral

La potencia es la probabilidad de rechazar H₀ cuando una alternativa específica es verdadera. Si β es la probabilidad de error tipo II para ese efecto, la potencia es 1 − β. No existe una única potencia independiente del efecto verdadero.

La «resolución» del ejemplo se formaliza como una diferencia que se desea poder detectar con una potencia y un nivel α determinados. Una diferencia mínima relevante debe justificarse por el fenómeno: no es la menor diferencia que una prueba pueda detectar alguna vez, ni una frontera que vuelva imposible detectar efectos menores.

Para planificar, σ representa una variabilidad supuesta o estimada a partir de información previa pertinente. Utilizar una estimación piloto introduce incertidumbre; conviene explorar varios escenarios. Los valores convencionales de tamaño del efecto no sustituyen una justificación sustantiva de δ.

Una muestra: detectar tres unidades

Se conservan δ = 3, σ = 10,2, potencia objetivo 0,90 y una alternativa bilateral. Aquí α = 0,05 se escribe explícitamente. El tamaño del efecto es d = 3/10,2 ≈ 0,29412. pwr es un paquete adicional, no una base de datos nativa de R.

# Instalacion opcional, una sola vez:
# install.packages("pwr")
# Si esta instalado, puede adjuntarse con library(pwr).

if (requireNamespace("pwr", quietly = TRUE)) {
  plan_pwr <- pwr::pwr.t.test(
    n = NULL, d = 3 / 10.2, sig.level = 0.05,
    power = 0.90, type = "one.sample", alternative = "two.sided"
  )
  print(plan_pwr)
  print(ceiling(plan_pwr$n))
}

El resultado es n ≈ 123,403, por lo que se necesitan 124 observaciones completas bajo el modelo y los supuestos adoptados. Con n = 124, la potencia calculada es aproximadamente 0,90139. Redondear al entero más cercano podría dejar la potencia por debajo del objetivo; se utiliza ceiling().

La función power.t.test() de R base permite resolver el mismo ejercicio sin instalar paquetes. Su argumento delta utiliza unidades originales y sd proporciona la desviación estándar. Se usa strict = TRUE para incluir, en el contraste bilateral, ambas regiones de rechazo y coincidir con el cálculo bilateral de pwr.^[12]

delta <- 3; sigma_plan <- 10.2; potencia_objetivo <- 0.9
d_plan <- delta / sigma_plan
plan_una <- power.t.test(n = NULL, delta = delta, sd = sigma_plan,
                        sig.level = alpha, power = potencia_objetivo,
                        type = "one.sample", alternative = "two.sided", strict = TRUE)
n_requerido <- ceiling(plan_una$n)
potencia_redondeada <- power.t.test(n = n_requerido, delta = delta,
                                  sd = sigma_plan, sig.level = alpha,
                                  type = "one.sample", strict = TRUE)$power
print(c(n = plan_una$n, entero = n_requerido,
        potencia = potencia_redondeada))

Con n = 40, ¿qué diferencia corresponde a potencia 0,90?

Al fijar n = 40 y resolver para δ, con la misma σ y α, resulta δ ≈ 5,3615 unidades. En la parametrización de pwr se resuelve primero para d y luego se multiplica por σ. «Óptimo» requiere una función de costos o un criterio adicional; aquí se está resolviendo una restricción de potencia, no un problema general de optimización.

plan_40 <- power.t.test(n = 40, delta = NULL, sd = sigma_plan,
                      sig.level = alpha, power = potencia_objetivo,
                      type = "one.sample", strict = TRUE)
print(plan_40$delta)
# Equivalente con pwr: resolver d = NULL y multiplicar por sigma_plan.

Una curva de planificación

Para cada diferencia entre dos y seis unidades se resuelve un tamaño muestral. Con α, σ y potencia fijos, detectar una diferencia menor requiere más observaciones. La curva puede presentarse continua antes del redondeo; la implementación final utiliza tamaños enteros y tiene forma escalonada.

deltas <- seq(2, 6, length.out = 30)
n_plan <- numeric(length(deltas)); d_grid <- deltas / sigma_plan
for (i in seq_along(deltas)) {
  n_plan[i] <- power.t.test(delta = deltas[i], sd = sigma_plan,
                          sig.level = alpha, power = potencia_objetivo,
                          type = "one.sample", strict = TRUE)$n
}
plot(deltas, ceiling(n_plan), type = "s", col = cobre, lwd = 2,
     xlab = "Diferencia a detectar", ylab = "Tamano muestral (redondeado)",
     main = "Potencia del 90%; alpha = 0.05")

Amilasa: dos grupos independientes

El segundo ejemplo conserva una varianza combinada de referencia de 290,9 y una diferencia de interés de 15 unidades/ml. La desviación estándar es √290,9 ≈ 17,056 y d ≈ 0,87947. Con α = 0,05 y potencia 0,90, el cálculo bilateral para dos grupos de igual tamaño da n ≈ 28,166 por grupo: se redondea a 29 por grupo, 58 en total.

s2p <- 290.9; sp <- sqrt(s2p); delta_amilasa <- 15
d_amilasa <- delta_amilasa / sp
plan_amilasa <- power.t.test(delta = delta_amilasa, sd = sp,
                           sig.level = alpha, power = potencia_objetivo,
                           type = "two.sample", strict = TRUE)
print(c(por_grupo = ceiling(plan_amilasa$n),
        total = 2 * ceiling(plan_amilasa$n)))

if (requireNamespace("pwr", quietly = TRUE)) {
  print(pwr::pwr.t.test(d = d_amilasa, sig.level = alpha,
                       power = 0.90, type = "two.sample"))
}

Este cálculo corresponde al modelo de dos grupos independientes con varianza común; no es un cálculo general de potencia de Welch con varianzas y tamaños arbitrarios. En muestras pareadas, n cuenta pares y la desviación estándar apropiada es la de las diferencias. Deben añadirse previsiones justificadas de pérdidas o registros incompletos; si se espera una fracción de retención r, una aproximación de reclutamiento es techo(n completo/r), siempre atendiendo al diseño.

Cómo se calcula β: la t no central

En la prueba de una media bajo normalidad, con efecto verdadero δ = μ − μ₀ y σ conocida como supuesto de planificación, T sigue una t no central con ν = n − 1 y no centralidad λ = δ√n/σ. La prueba sigue estimando el error estándar con S; conocer σ como supuesto de potencia no convierte el contraste en una prueba z.

gl_plan <- n_requerido - 1
crit_plan <- qt(1 - alpha/2, gl_plan)
ncp_plan <- delta * sqrt(n_requerido) / sigma_plan
beta_plan <- pt(crit_plan, gl_plan, ncp = ncp_plan) -
             pt(-crit_plan, gl_plan, ncp = ncp_plan)
print(c(beta = beta_plan, potencia = 1 - beta_plan))
stopifnot(abs(1 - beta_plan - potencia_redondeada) < 1e-10)

Con n = 124 y los parámetros anteriores se obtiene β ≈ 0,09861 y potencia ≈ 0,90139. La concordancia con la función de planificación permite comprobar el fundamento del cálculo. El script verifica esa igualdad numéricamente.

13 / Criterio científico

Supuestos, límites y comunicación de resultados

Una función devuelve un número aun cuando la pregunta esté mal formulada. La calidad del análisis depende de la relación entre diseño, modelo, medición y fenómeno. Los siguientes criterios organizan esa revisión sin convertirla en una secuencia mecánica de pruebas preliminares.

Un gráfico cuantil–cuantil y la inspección de los datos pueden ayudar a evaluar la forma de la distribución. No rechazar una prueba de normalidad no demuestra normalidad, especialmente con muestras pequeñas. Tampoco existe una obligación de rechazar todos los métodos paramétricos por una desviación menor: importa la sensibilidad de la inferencia concreta a esa desviación.

Si se cambian las herramientas, debe conservarse la pregunta o declarar su modificación. Una prueba de rangos no es automáticamente una prueba de igualdad de medias. Un remuestreo bootstrap o una permutación requieren condiciones propias; permutar etiquetas exige una justificación de intercambiabilidad y no repara por sí solo un diseño dependiente o sesgado.

Significación, geometría y magnitud

Es posible representar geométricamente las observaciones como vectores y las medias como proyecciones sobre la dirección constante. En ese lenguaje, una diferencia entre estimaciones tiene magnitud y dirección. Sin embargo, la significación estadística requiere además una distribución de referencia, una regla de extremidad y supuestos de muestreo. No se define únicamente como una variación visible del patrón geométrico.

Un efecto pequeño puede resultar significativo con gran precisión; uno importante puede quedar acompañado de un intervalo amplio. La decisión binaria no agota la evidencia. Si se realizan numerosos contrastes, hay que definir la familia de preguntas y considerar el control de multiplicidad apropiado; por ejemplo, p.adjust(valores_p, method = "holm") para un ajuste de Holm de una familia determinada. No se necesita ajustar indiscriminadamente todo cálculo presentado con finalidad didáctica.

Una forma de informar el resultado

Un informe debe indicar pregunta y población, diseño, tamaños realmente analizados, estimación en unidades originales, intervalo, método y supuestos, estadístico, grados de libertad y valor p, junto a las limitaciones. Para el ejemplo de sueño: «La diferencia media A − B fue −1,58 horas (IC del 95 %: −2,460 a −0,700); prueba t pareada, t(9) = −4,062, p = 0,002833, n = 10 pares». Después debe explicarse qué significa esa magnitud y hasta dónde permite generalizar el diseño.

Un valor p muy pequeño se informa con cifras suficientes o notación científica, no como p = 0. Un resultado no significativo se acompaña de su estimación e incertidumbre; no se reformula como inexistencia del efecto. Una afirmación de equivalencia exige márgenes y un análisis específico.

14 / Consulta

Guía de sintaxis y documentos reproducibles

La inspección de un objeto de prueba con names(prueba_una) permite localizar sus componentes. prueba_una$p.value, $statistic, $parameter, $conf.int y $estimate extraen resultados sin copiarlos manualmente de la consola. El nombre y la interpretación de cada componente deben verificarse en la documentación de la función.

R Markdown: resultados dentro de la prosa

En un archivo .Rmd, una expresión R entre comillas invertidas y precedida por r se evalúa al compilar el documento. No es una expresión que WordPress ejecute ni una instrucción que deba pegarse tal cual en la consola.

La media de incremento de sueno es `r mean(sleep$extra)` horas.
Hay `r nrow(sleep)` registros y `r length(unique(sleep$ID))` pacientes.

Así se distinguen veinte registros y diez pacientes en sleep. La media de sleep$extra mezcla los dos tratamientos; para una comparación se necesitan sus medias por tratamiento o, en este diseño, las diferencias por paciente. Un documento reproducible debe registrar también las transformaciones, los criterios de exclusión y las versiones del software.

15 / Apéndice ejecutable

Script completo reproducible

El siguiente script reúne todos los pasos en orden, incluye los 21 registros ilustrativos, produce los gráficos y compara los cálculos manuales con las funciones de R mediante comprobaciones numéricas. Se ejecuta en una sesión limpia. Los apartados principales funcionan con R base; si pwr y gplots están instalados, se ejecutan también las variantes opcionales. No se instalan paquetes automáticamente.

Para utilizar el archivo original cuando se disponga de él, cambie usar_archivo_original a TRUE y adapte su importación al formato real. Los resultados impresos de tiendas cambiarán; los del conjunto histórico sleep y los escenarios de planificación no dependen de esos registros.

Descargar el script en R

# FUNDAMENTOS DE ESTIMACION Y PRUEBAS DE HIPOTESIS EN R
# JOSÉ MAURICIO GÓMEZ JULIÁN
# Edicion revisada. Los 21 registros son ILUSTRATIVOS, no datos recuperados.
# Ejecutar en una sesion limpia. No se instalan paquetes automaticamente.

# 1. Datos y comprobaciones
usar_archivo_original <- FALSE
if (usar_archivo_original) {
  estudios <- read.table("estudios.txt", header = FALSE,
                        col.names = c("menor16", "ipc", "ventas"))
} else {
  estudios <- data.frame(
    menor16 = c(12,15,11,18,14,16,13,19,17,12,20,15,18,14,16,19,13,17,11,20,18),
    ipc = c(11,12,14,16,13,15,10,16,14,12,18,20,17,22,19,24,18,21,23,17,20),
    ventas = c(145,158,172,183,167,192,151,176,164,180,188,205,194,217,201,229,182,212,221,196,208)
  )
}
stopifnot(all(c("menor16", "ipc", "ventas") %in% names(estudios)),
          all(vapply(estudios, is.numeric, logical(1))))
print(names(estudios)); print(nrow(estudios)); print(ncol(estudios))
print(dim(estudios)); str(estudios)
print(colSums(is.na(estudios)))
if (any(!is.finite(as.matrix(estudios)))) {
  stop("Hay valores ausentes o no finitos: documente su tratamiento antes de analizar.")
}
stopifnot(nrow(estudios) >= 2L)
ventas <- estudios$ventas
print(c(media = mean(ventas), sd = sd(ventas), varianza = var(ventas)))
print(apply(estudios, 2, mean)); print(apply(estudios, 2, sd))
print(quantile(ventas)); print(apply(estudios, 2, quantile))
qv <- quantile(ventas, probs = c(0.025, 0.975))
print(qv)

# Colores y graficos base: no requieren paquetes adicionales.
tinta <- "#2c2418"; cobre <- "#8b4513"; beige <- "#d4c8b0"
par(col = tinta, col.axis = tinta, col.lab = tinta, col.main = tinta,
    bg = "#fffdf7", las = 1)
hist(ventas, breaks = 7, col = beige, border = cobre,
     xlab = "Ventas", ylab = "Frecuencia", main = "Ventas: datos ilustrativos")
abline(v = qv, col = cobre, lty = 2)

# 2. Simulacion: parametros poblacionales, cuantiles y densidad
set.seed(20210906)
y <- rnorm(1000, mean = 35, sd = 2)
qy <- quantile(y, probs = c(0.025, 0.975))
hist(y, probability = TRUE, col = beige, border = cobre,
     main = "N(35, 4): densidad estimada y teorica", xlab = "Valor")
abline(v = qy, col = cobre, lty = 2)
lines(density(y), col = cobre, lwd = 2)
curve(dnorm(x, mean = 35, sd = 2), add = TRUE, lty = 3, lwd = 2)
legend("topright", c("Estimada por nucleos", "Normal teorica"),
       col = c(cobre, tinta), lty = c(1, 3), bty = "n", cex = 0.7)

plot.m <- function(n, m, s) {
  stopifnot(length(n) == 1L, is.finite(n), n >= 2, n == as.integer(n),
            length(m) == 1L, is.finite(m), length(s) == 1L, is.finite(s), s > 0)
  y <- rnorm(n, mean = m, sd = s)
  qy <- quantile(y, c(0.025, 0.975))
  hist(y, probability = TRUE, col = beige, border = cobre,
       main = paste("Muestra normal: n =", n), xlab = "Valor")
  abline(v = qy, col = cobre, lty = 2)
  abline(v = mean(y), col = tinta, lwd = 2)
  lines(density(y), col = cobre)
  invisible(list(media = mean(y), sd = sd(y), cuantiles = qy, datos = y))
}
muestra_10 <- plot.m(10, 100, 15)
muestra_10000 <- plot.m(10000, 100, 15)
print(muestra_10[c("media", "sd", "cuantiles")])
print(muestra_10000[c("media", "sd", "cuantiles")])

# 3. Distribucion muestral de la media: bucle y equivalencia teorica
n_sim <- 10000L; mu_sim <- 100; sigma_sim <- 15; B <- 1000L
medias <- numeric(B)
for (i in seq_len(B)) {
  sam <- rnorm(n_sim, mean = mu_sim, sd = sigma_sim)
  medias[i] <- mean(sam)
}
print(c(media_teorica = mu_sim, media_simulada = mean(medias),
        ee_teorico = sigma_sim / sqrt(n_sim), ee_simulado = sd(medias)))
qm <- quantile(medias, c(0.025, 0.975))
hist(medias, probability = TRUE, col = beige, border = cobre,
     main = "Distribucion de medias", xlab = "Media muestral")
abline(v = qm, col = cobre, lty = 2); lines(density(medias), col = cobre)

# Extension: poblacion exponencial de media 1 y desviacion estandar 1.
set.seed(20210907)
medias_exp <- lapply(c(1L, 5L, 30L), function(nn) {
  replicate(4000L, mean(rexp(nn, rate = 1)))
})
names(medias_exp) <- c("n1", "n5", "n30")
print(vapply(medias_exp, mean, numeric(1)))
print(vapply(medias_exp, sd, numeric(1)))

# 4. Intervalos de confianza: funcion con validaciones
ic <- function(x, p = 0.95, na.rm = FALSE) {
  stopifnot(is.numeric(x), length(p) == 1L, is.finite(p), p > 0, p < 1,
            is.logical(na.rm), length(na.rm) == 1L, !is.na(na.rm))
  if (na.rm) x <- x[!is.na(x)]
  if (any(!is.finite(x))) stop("Datos ausentes o no finitos.")
  n_x <- length(x)
  if (n_x < 2L) stop("Se requieren al menos dos observaciones.")
  if (sd(x) == 0) stop("No hay variacion muestral.")
  ee_x <- sd(x) / sqrt(n_x)
  critico <- qt((1 + p) / 2, df = n_x - 1)
  setNames(mean(x) + c(-1, 1) * critico * ee_x, c("inferior", "superior"))
}
print(ic(ventas, 0.95)); print(ic(ventas, 0.99))
print(ic(estudios$ipc, 0.95)); print(ic(estudios$menor16, 0.95))
stopifnot(isTRUE(all.equal(unname(ic(ventas)),
                          as.numeric(t.test(ventas)$conf.int))))

# 5. Prueba t de una muestra: calculo manual y automatico
alpha <- 0.05; mu0 <- 180
n_ventas <- length(ventas); gl <- n_ventas - 1L
ee_ventas <- sd(ventas) / sqrt(n_ventas)
t_obs <- (mean(ventas) - mu0) / ee_ventas
p_manual <- 2 * pt(-abs(t_obs), df = gl)
prueba_una <- t.test(ventas, mu = mu0, alternative = "two.sided",
                    conf.level = 1 - alpha)
print(prueba_una)
print(c(t_manual = t_obs, p_manual = p_manual))
stopifnot(isTRUE(all.equal(p_manual, prueba_una$p.value)))
# Una alternativa unilateral se elige antes de mirar los resultados.
prueba_mayor <- t.test(ventas, mu = mu0, alternative = "greater")

# 6. Agrupacion y resumen por grupo. Umbral ilustrativo, no optimizado.
umbral <- 17
estudios$ipc1 <- ifelse(estudios$ipc < umbral, 1L, 2L)
estudios$ipc2 <- factor(estudios$ipc1, levels = c(1, 2),
                        labels = c("ipc < 17", "ipc >= 17"))
plot(ipc ~ ipc2, data = estudios, col = beige, border = cobre,
     xlab = "Grupo", ylab = "IPC", main = "Clasificacion por IPC")
abline(h = umbral, col = cobre, lty = 2)
boxplot(ventas ~ ipc2, data = estudios, col = beige, border = cobre,
        ylab = "Ventas", xlab = "Grupo", main = "Ventas por grupo")
med_grupos <- tapply(estudios$ventas, estudios$ipc2, mean)
sd_grupos <- tapply(estudios$ventas, estudios$ipc2, sd)
var_grupos <- tapply(estudios$ventas, estudios$ipc2, var)
n_grupos <- table(estudios$ipc2)
stopifnot(length(n_grupos) == 2L, all(n_grupos >= 2L), all(var_grupos > 0))
print(cbind(n = as.numeric(n_grupos), media = med_grupos,
            sd = sd_grupos, varianza = var_grupos))
ic_grupos <- sapply(split(estudios$ventas, estudios$ipc2), ic)
plot(seq_along(med_grupos), med_grupos, pch = 18, xaxt = "n",
     xlim = c(0.5, 2.5), ylim = range(ic_grupos), xlab = "Grupo", ylab = "Ventas",
     main = "Medias e intervalos de confianza del 95%")
axis(1, at = 1:2, labels = names(med_grupos))
arrows(1:2, ic_grupos[1, ], 1:2, ic_grupos[2, ],
       angle = 90, code = 3, length = 0.05, col = cobre)
# Recurso opcional del articulo original:
if (requireNamespace("gplots", quietly = TRUE)) {
  gplots::plotmeans(ventas ~ ipc2, data = estudios, connect = FALSE,
                    col = cobre, barcol = cobre)
}

# 7. Medias independientes: varianza combinada y Welch
dif_medias <- unname(med_grupos[1] - med_grupos[2])
sp2 <- sum((as.numeric(n_grupos) - 1) * var_grupos) / (sum(n_grupos) - 2)
ee_pool <- sqrt(sp2 * sum(1 / as.numeric(n_grupos)))
t_pool <- dif_medias / ee_pool
gl_pool <- sum(n_grupos) - 2
p_pool <- 2 * pt(-abs(t_pool), df = gl_pool)
prueba_pool <- t.test(ventas ~ ipc2, data = estudios, var.equal = TRUE)
a_w <- var_grupos / as.numeric(n_grupos)
ee_welch <- sqrt(sum(a_w))
gl_welch <- sum(a_w)^2 / sum(a_w^2 / (as.numeric(n_grupos) - 1))
t_welch <- dif_medias / ee_welch
p_welch <- 2 * pt(-abs(t_welch), df = gl_welch)
prueba_welch <- t.test(ventas ~ ipc2, data = estudios, var.equal = FALSE)
print(prueba_pool); print(prueba_welch)
stopifnot(isTRUE(all.equal(unname(p_pool), prueba_pool$p.value)),
          isTRUE(all.equal(unname(p_welch), prueba_welch$p.value)),
          isTRUE(all.equal(unname(t.test(ventas ~ ipc2, data = estudios)$p.value),
                           prueba_welch$p.value)))

# 8. Prueba F e intervalo de la razon de varianzas
razon0 <- 1
razon_muestral <- unname(var_grupos[1] / var_grupos[2])
f_obs <- razon_muestral / razon0
gl1 <- as.numeric(n_grupos[1]) - 1; gl2 <- as.numeric(n_grupos[2]) - 1
p_f <- min(1, 2 * min(pf(f_obs, gl1, gl2),
                       pf(f_obs, gl1, gl2, lower.tail = FALSE)))
ic_razon <- razon_muestral / qf(c(1 - alpha/2, alpha/2), gl1, gl2)
prueba_f <- var.test(ventas ~ ipc2, data = estudios, ratio = razon0)
print(prueba_f); print(ic_razon)
stopifnot(isTRUE(all.equal(p_f, prueba_f$p.value)),
          isTRUE(all.equal(ic_razon, as.numeric(prueba_f$conf.int))))

# 9. Muestras pareadas: alinear por ID, no por posicion supuesta
data("sleep", package = "datasets")
A <- sleep[sleep$group == "1", c("ID", "extra")]
B_sleep <- sleep[sleep$group == "2", c("ID", "extra")]
stopifnot(!anyDuplicated(A$ID), !anyDuplicated(B_sleep$ID),
          setequal(A$ID, B_sleep$ID))
pares <- merge(A, B_sleep, by = "ID", suffixes = c("_A", "_B"), sort = TRUE)
stopifnot(nrow(pares) == 10L, all(complete.cases(pares)))
dif_sueno <- pares$extra_A - pares$extra_B
prueba_pareada <- t.test(pares$extra_A, pares$extra_B, paired = TRUE)
prueba_diferencias <- t.test(dif_sueno, mu = 0)
prueba_independiente <- t.test(pares$extra_A, pares$extra_B)
print(prueba_pareada); print(prueba_independiente)
print(c(dz = mean(dif_sueno) / sd(dif_sueno)))
stopifnot(isTRUE(all.equal(prueba_pareada$p.value, prueba_diferencias$p.value)))
matplot(1:2, t(as.matrix(pares[c("extra_A", "extra_B")])),
        type = "b", pch = c(1, 16), lty = 1, col = cobre,
        xaxt = "n", xlab = "Tratamiento", ylab = "Incremento de horas de sueno",
        main = "Mediciones pareadas por paciente")
axis(1, at = 1:2, labels = c("A", "B"))
# Grafico por paciente conservado del ejemplo original:
plot(seq_len(nrow(pares)), pares$extra_A, type = "b", pch = 1,
     col = cobre, ylim = range(pares[c("extra_A", "extra_B")]),
     xlab = "Paciente", ylab = "Incremento de horas", main = "A y B por paciente")
lines(seq_len(nrow(pares)), pares$extra_B, type = "b", pch = 16, lty = 2)
legend("topleft", c("A", "B"), col = c(cobre, tinta), pch = c(1, 16),
       lty = c(1, 2), bty = "n")

# 10. Potencia y tamano muestral, usando R base
delta <- 3; sigma_plan <- 10.2; potencia_objetivo <- 0.9
d_plan <- delta / sigma_plan
plan_una <- power.t.test(n = NULL, delta = delta, sd = sigma_plan,
                        sig.level = alpha, power = potencia_objetivo,
                        type = "one.sample", alternative = "two.sided", strict = TRUE)
n_requerido <- ceiling(plan_una$n)
potencia_redondeada <- power.t.test(n = n_requerido, delta = delta,
                                  sd = sigma_plan, sig.level = alpha,
                                  type = "one.sample", strict = TRUE)$power
plan_40 <- power.t.test(n = 40, delta = NULL, sd = sigma_plan,
                      sig.level = alpha, power = potencia_objetivo,
                      type = "one.sample", strict = TRUE)
deltas <- seq(2, 6, length.out = 30)
n_plan <- numeric(length(deltas)); d_grid <- deltas / sigma_plan
for (i in seq_along(deltas)) {
  n_plan[i] <- power.t.test(delta = deltas[i], sd = sigma_plan,
                          sig.level = alpha, power = potencia_objetivo,
                          type = "one.sample", strict = TRUE)$n
}
plot(deltas, ceiling(n_plan), type = "s", col = cobre, lwd = 2,
     xlab = "Diferencia a detectar", ylab = "Tamano muestral (redondeado)",
     main = "Potencia del 90%; alpha = 0.05")
# Amilasa: dos grupos independientes, iguales tamanos y varianza comun.
s2p <- 290.9; sp <- sqrt(s2p); delta_amilasa <- 15
d_amilasa <- delta_amilasa / sp
plan_amilasa <- power.t.test(delta = delta_amilasa, sd = sp,
                           sig.level = alpha, power = potencia_objetivo,
                           type = "two.sample", strict = TRUE)
print(c(d = d_plan, n_continuo = plan_una$n, n_requerido = n_requerido,
        potencia_redondeada = potencia_redondeada,
        delta_con_n40 = plan_40$delta, d_amilasa = d_amilasa,
        n_amilasa_por_grupo = plan_amilasa$n,
        n_amilasa_redondeado = ceiling(plan_amilasa$n)))

# Version original con paquete pwr, opcional:
# install.packages(c("pwr", "gplots")) # ejecutar solo si se desea instalarlos
if (requireNamespace("pwr", quietly = TRUE)) {
  plan_pwr <- pwr::pwr.t.test(n = NULL, d = d_plan, sig.level = alpha,
                            power = potencia_objetivo, type = "one.sample",
                            alternative = "two.sided")
  plan_pwr40 <- pwr::pwr.t.test(n = 40, d = NULL, sig.level = alpha,
                              power = potencia_objetivo, type = "one.sample")
  plan_pwr_amilasa <- pwr::pwr.t.test(n = NULL, d = d_amilasa,
                                   sig.level = alpha, power = potencia_objetivo,
                                   type = "two.sample")
  print(plan_pwr); print(plan_pwr40$d * sigma_plan); print(plan_pwr_amilasa)
  stopifnot(abs(plan_pwr$n - plan_una$n) < 0.001,
            abs(plan_pwr_amilasa$n - plan_amilasa$n) < 0.001)
}

# 11. Potencia exacta para una media bajo normalidad: t no central
gl_plan <- n_requerido - 1
crit_plan <- qt(1 - alpha/2, gl_plan)
ncp_plan <- delta * sqrt(n_requerido) / sigma_plan
beta_plan <- pt(crit_plan, gl_plan, ncp = ncp_plan) -
             pt(-crit_plan, gl_plan, ncp = ncp_plan)
print(c(beta = beta_plan, potencia = 1 - beta_plan))
stopifnot(abs(1 - beta_plan - potencia_redondeada) < 1e-10)

# R Markdown: esto se escribe dentro de un .Rmd, no como codigo R aislado:
# La media de incremento de sueno es `r mean(sleep$extra)` horas.
# El numero de registros es `r nrow(sleep)`; hay `r length(unique(sleep$ID))` pacientes.

# 12. Resumen verificable y version de la sesion
cat("\nRESUMEN_VERIFICADO\n")
print(list(
  media_ventas = mean(ventas), sd_ventas = sd(ventas), cuantiles_ventas = qv,
  ic95_ventas = ic(ventas), ic99_ventas = ic(ventas, 0.99),
  una_muestra = c(t = unname(prueba_una$statistic), gl = gl, p = prueba_una$p.value),
  grupos = c(diferencia = dif_medias, t_welch = unname(t_welch),
             gl_welch = unname(gl_welch), p_welch = p_welch, p_pool = p_pool,
             f = f_obs, p_f = p_f),
  sueno = c(diferencia = mean(dif_sueno), sd_diferencias = sd(dif_sueno),
            t = unname(prueba_pareada$statistic), p = prueba_pareada$p.value,
            ic_inf = prueba_pareada$conf.int[1], ic_sup = prueba_pareada$conf.int[2],
            p_independiente = prueba_independiente$p.value),
  planificacion = c(n_una = plan_una$n, delta40 = plan_40$delta,
                    n_amilasa = plan_amilasa$n)))
sessionInfo()

16 / Fuentes

Referencias y documentación

La documentación enlazada permite comprobar argumentos, distribuciones y supuestos. Las funciones se citan donde se utilizan; los enlaces de consulta no sustituyen las explicaciones desarrolladas en el texto.

1. Wasserstein, R. L., y Lazar, N. A. (2016). The ASA’s Statement on p-Values: Context, Process, and Purpose. The American Statistician, 70(2), 129–133. doi:10.1080/00031305.2016.1154108.
2. R Core Team. Apply Functions Over Array Margins: apply().
3. R Core Team. Sample Quantiles: quantile().
4. R Core Team. The Normal Distribution: dnorm(), pnorm(), qnorm(), rnorm().
5. R Core Team. Kernel Density Estimation: density().
6. R Core Team. The Student t Distribution: dt(), pt(), qt(), rt().
7. R Core Team. Apply a Function Over a Ragged Array: tapply().
8. R Core Team. Student’s t-Test: t.test().
9. R Core Team. The F Distribution: df(), pf(), qf(), rf().
10. R Core Team. F Test to Compare Two Variances: var.test(). Véase también fligner.test().
11. R Core Team. Student’s Sleep Data: sleep. La documentación remite a Cushny y Peebles (1905) y Student (1908).
12. R Core Team. Power Calculations for One and Two Sample t Tests: power.t.test().
13. Champely, S., y colaboradores. pwr: Basic Functions for Power Analysis, versión 1.3-0. Manual del paquete. Los cálculos se basan en la parametrización de tamaños del efecto de Cohen (1988).
14. Equipo de gplots. gplots: Various R Programming Tools for Plotting Data. Manual de plotmeans() y otras funciones.
15. R Core Team. Histogram: hist(); Adjust P-values for Multiple Comparisons: p.adjust().

Publicación original y lecturas que acompañaban sus comentarios

Esta edición revisa Fundamentos generales del proceso de estimación y prueba de hipótesis en R Studio. Parte II, código en R Studio con comentarios (2021), publicada bajo el seudónimo Isadore Nabi. La firma autoral de esta edición es JOSÉ MAURICIO GÓMEZ JULIÁN.

Se conservan como referencias de procedencia los recursos citados en los comentarios originales: Minitab, explicación de valores t y pruebas t; Statology, distribución t en R; Boston University, potencia y tamaño muestral; y StatTrek, calculadora de la distribución F. Se incluyen como parte de la trayectoria del texto; la documentación oficial fundamenta las correcciones de sintaxis y funciones de esta edición.

También se mantienen las conexiones temáticas con los textos del autor Análisis teórico de la función cuantil en R Studio y Modelos lineales generalizados. La Parte I de teoría estadística es mencionada en la publicación original; este cuaderno desarrolla aquí los fundamentos utilizados sin atribuirle contenido que no se haya incorporado expresamente.