Errores tipo M y tipo S
jmm
4/10/2019
Nuestro héroe Andrew Gelman nunca no está preguntándose cuánto podemos realmente confiar en los resultados de nuestros análisis de datos.
Un ejemplo claro de este problema son los llamados errores tipo M y S. La idea es la siguiente. Si un efecto (real) es de poca magnitud y tenemos muestras chicas o ruidosas (mediciones con poca precisión), solamente vamos a detectar efectos estadísticamente significativos cuando la magnitud estimada del efecto es exagerada (error M) y esto puede ocurrir incluso con efectos de signo contrario al verdadero (error S). Una buena referencia para estos errores se puede ver aquí. O buscando en el blog de Gelman https://statmodeling.stat.columbia.edu/
Veamos de qué se trata haciendo unas simulaciones. Vamos a suponer que existe un efecto real entre una variable predictora y una respuesta, pero que la magnitud del efecto no es muy grande. Podemos simlular el caso de una regresión simple:
set.seed(123)
n = 20
b0 = 0.5
b1 = 0.1
sigma = 1
x = runif(n, -3,3)
y = rnorm(n, b0 + b1 * x, sd = sigma)
plot(x, y)
curve(b0 + b1 * x, add = TRUE, lwd = 2, col = 2)
abline(lm(y~x))
Ejemplo de datos simulados. La línea roja muestra la relación verdadera entre la variable predictora y la media de la variable de respuesta. La línea negra muestra la regresión para el set de datos simulados.
Podemos hacer unas cuantas réplicas y ver cuándo tenemos estimaciones significativas para la pendiente (b1) y luego graficar las relaciones estimadas entre \(x\) e \(y\)
n.reps = 1000
n = 20
res = matrix(NA, n.reps, 2)
se = numeric(n.reps) # para guardar los errores estándar de la pendiente
ps = numeric(n.reps) # para guardar los valores de p de la pendiente
for(i in 1: n.reps){
x = runif(n, -3,3)
y = rnorm(n, b0 + b1 * x, sd = sigma)
tmp = summary(lm(y~x))$coefficients
se[i] = tmp[2,2]
ps[i] = tmp[2,4]
res[i, ] = tmp[,1]
}Grafiquemos las regresiones que fueron significativas
curve(b0 + b1 * x, xlim = c(-3,3), ylim = c(-0.5, 3))
for(i in 1: n.reps){
if(ps[i] <= 0.05){
curve(res[i,1] + res[i,2] * x, add = TRUE, col=rgb(0,0,0, 0.2))
}
}
curve(b0 + b1 * x, add = TRUE, col = 2, lwd = 3)
Errores tipo M y S. Las líneas grises muestran las relaciones estimadas entre la variable predictora y el valor esperado de la variable de respuesta cuando la pendiente resultó significativa. En rojo se muestra la relación verdadera. En general, si la pendiente es significativa, la magnitu del efeto está exagerada (error M) e incluso puede ser del signo contrario al verdadero (error S)
Para calcular la magnitud de la tasa de exageración (error M) hacemos:
tmp <- which(ps <= 0.05)
mean(abs(res[tmp,2])/b1)## [1] 3.136809Quiere decir que, en promedio, cuando la pendiente es significativa, los tamaños de efecto estimados serán al menos tres veces mayores al verdadero.
También podemos calcular el poder (probabilidad de rechazar correctamente la hipótesis nula) de este análisis a partir de el valor de significancia buscado, los grados de libertad, la magnitud del efecto verdadero y el error estándar de la estimación:
alpha = 0.05
A = 0.1
df = 498
s = mean(se)
z <- qt(1-alpha/2, df)
p.hi <- 1 - pt(z-A/s, df)
p.lo <- pt(-z-A/s, df)
power <- p.hi + p.lo
power## [1] 0.1162243Ejercicios
- Cuántas observaciones necesito para obtener un error M de menos de 1.2?
- Cuánto tengo que reducir sigma para obtener un error M de menos de 1.2 con un tamaño muestral de 20?
Análisis Bayesianos y rol de las previas como regularizadoras
Ahora podemos ver qué pasa si hacemos análisis bayesianos con diferentes previas. Para estos análisis vamos a escribir el modelo de regresión en Stan:
data {
int<lower=1> N;
vector[N] Y;
vector[N] X;
real<lower=0> sd_b; // sd prior for slope
}
parameters {
real b0;
real b1;
real<lower=0> sigma;
}
model {
vector[N] mu = b0 + b1 * X;
b0 ~ normal(0,1);
b1 ~ normal(0, sd_b);
sigma ~ student_t(3,0,10);
Y ~ normal(mu, sigma);
}
Hagamos una regresión para ver cómo funciona:
library(rstan)## Loading required package: ggplot2## Loading required package: StanHeaders## rstan (Version 2.18.2, GitRev: 2e1f913d3ca3)## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
stan.data = list(N = n,
Y = y,
X = x,
sd_b = 0.2)
fit <- stan(file = "reg.stan",
data = stan.data,
iter = 1000,
chains = 3)
print(fit)## Inference for Stan model: reg.
## 3 chains, each with iter=1000; warmup=500; thin=1;
## post-warmup draws per chain=500, total post-warmup draws=1500.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## b0 0.59 0.01 0.21 0.17 0.45 0.59 0.73 1.01 1126 1
## b1 0.14 0.00 0.11 -0.09 0.07 0.14 0.22 0.35 1367 1
## sigma 0.96 0.01 0.18 0.68 0.83 0.93 1.06 1.37 746 1
## lp__ -8.96 0.06 1.33 -12.36 -9.56 -8.60 -8.01 -7.45 522 1
##
## Samples were drawn using NUTS(diag_e) at Wed Jun 12 21:39:35 2019.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).Ahora deberíamos hacer varias réplicas para ver en qué casos obtenemos errores tipo M y S según las previas que usamos. Para ahorrar tiempo vamos a aproximar las posteriores usando la funcion map (maximum a posteriori) del paquete rethinking
library(rethinking)## Loading required package: parallel## rethinking (Version 1.88)m <- map(
alist(
y ~ dnorm( mu , sigma ) ,
mu <- b0 + b1*x ,
b0 ~ dnorm( 0 , 1 ) ,
b1 ~ dnorm( 0 , 0.5 ) ,
sigma ~ dunif( 0 , 5 )
), data = data.frame(x,y) )n.reps = 1000
resB = matrix(NA, n.reps, 2)
seB = numeric(n.reps) # para guardar los errores estándar de la pendiente
psB = numeric(n.reps) # posteriores de b1 distintas de cero?
for(i in 1: n.reps){
x = runif(n, -3,3)
y = rnorm(n, b0 + b1 * x, sd = sigma)
start <- list(
b0=mean(y),
b1=0,
sigma = 1
)
m <- map(
alist(
y ~ dnorm( mu , sigma ) ,
mu <- b0 + b1*x ,
b0 ~ dnorm( 0 , 1 ) ,
b1 ~ dnorm( 0 , 0.25 ) ,
sigma ~ dunif( 0 , 5 )
),
data = data.frame(x,y),
start = start)
tmp = precis(m, prob = 0.95)
resB[i, ] = tmp[c("b0", "b1"),"mean"]
seB[i] = tmp["b1", "sd"]
if(sign(tmp["b1", "2.5%"])==sign(tmp["b1", "97.5%"])){
psB[i] = 1
}
}Hacemos la gráfica
curve(b0 + b1 * x, xlim = c(-3,3), ylim = c(-0.5, 3))
for(i in 1: n.reps){
if(psB[i] == 1){
curve(resB[i,1] + resB[i,2] * x, add = TRUE, col=rgb(0,0,0, 0.2))
}
}
curve(b0 + b1 * x, add = TRUE, col = 2, lwd = 3)
Y calculamos la magnitud de la tasa de exageración (error M) como
tmp <- which(psB==1)
mean(abs(resB[tmp,2])/b1)## [1] 2.450643Vemos que al usar una previa que restringe la posible magnitud de la pendiente, reducimos bastante la tasa de exageración.
Ejercicios
Armarse de paciencia y repetir el ejemplo de arriba pero usando stan en lugar de map. Para eso pueden usar el siguiente script
library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
n.reps = 1000
resB = matrix(NA, n.reps, 2)
seB = numeric(n.reps)
psB = numeric(n.reps)
for(i in 1: n.reps){
x = runif(n, -3,3)
y = rnorm(n, b0 + b1 * x, sd = sigma)
stan.data = list(N = n,
Y = y,
X = x,
sd_b = 1)
fit <- stan(file = "reg.stan",
data = stan.data,
iter = 1000,
chains = 3,
refresh = 0)
tmp = summary(fit)$summary
if(all(tmp[,10] < 1.1)){
resB[i, ] = tmp[1:2,1]
seB[i] = tmp[2,3]
if(sign(tmp[2,4])==sign(tmp[2,8])){
psB[i] = 1
}
}
}- Comparen la tasa de exageración para
sd_bde \(1\), \(0.2\) y \(0.1\)