Estudio de asociación genética con un locus multialélico: Estimación del poder estadístico mediante simulaciones
Diego Véliz
23 de abril de 2017
Se compara el poder estadístico de 4 pruebas diferentes: La prueba Chi-cuadrado general, que evalúa la independencia entre la enfermedad y los genotipos. En esta prueba, los genotipos se clasificaron en dos categorías: Presenta el alelo E4 vs no presenta el alelo E4.
Además, se evalúa la prueba de Cochran-Armitage Trend Test (CATT), empleada para datos ordinales (el efecto de los genotipos va en el orden AA>Aa>aa). Se ordenaron los genotipos de acuerdo al número de copias del alelo E4, es decir: E4E4 > E3E4/E2E4 > E2E2/E2E3/E3E3.
También se incluye un modelo de regresión logística empleando el número de copias del alelo E4 como variable predictora.
Finalmente, se incluye la prueba de asociación alélica, que compara el número de copias del alelo E4 en población afectada vs población control. Esta prueba es particularmente potente cuando el RR de la enfermedad sigue un modelo multiplicativo, pero requiere HWE.
Notar que la frecuencia de los genotipos E4E4 es muy baja en comparación con la frecuencia del genotipo E3E4, por lo tanto, el poder estadístico depende principalmente del RR del genotipo E3E4.
En la siguiente tabla se presenta la equivalencia entre los valores de RR y OR para los parámetros empleados en esta simulación.
Los datos simulados fueron generados con los RR de esta tabla, pero el modelo de asociación empleado sobre los datos generados fueron los mencionados anteriormente.
EQUIVALENCIA ENTRE LOS VALORES EMPLEADOS DE RIESGO RELATIVO Y ODDS RATIO ESPERADO
| Relative Risk | Odds Ratio | |
|---|---|---|
| E2E2/E2E3 | 0.8 | 0.80 |
| E3E3 | 1.0 | 1.00 |
| E2E4/E3E4 | 3.0 | 3.16 |
| E4/E4 | 5.0 | 5.58 |
| E4/E4 | 7.0 | 8.27 |
| E4/E4 | 8.0 | 9.72 |
| E4/E4 | 9.0 | 11.27 |
| E4/E4 | 10.0 | 12.90 |
PARÁMETROS EMPLEADOS EN LA SIMULACIÓN
Tamaños muestrales: 50, 100, 150, 200, 250, 300
prevalencia = 0.0345
Proporción de casos con respecto al total de las muestras: 0.5
Número de réplicas empleadas en la simulación: 10000 (a mayor número de réplicas, mayor exactitud). Para 1000 réplicas la máxima desviación estandar del poder estadístico estimado es 0.0158 (1.58%). Para 10000 réplicas la desviación estándar máxima es 0.005 (0.5%)
Valores de Riesgo Relativo del genotipo E4E4: 5, 7, 8, 9 y 10
Método de análisis: Se comparan 4 modelos diferentes: Chi-cuadrado, Cochran-Armitage Trend Test, Regresión logística y asociación alélica. Los dos últimos métodos se detallan en el libro “Statistical Genetics: Gene mapping through linkage and associaction.”
Frecuencia del alelo E2: 0.1
Frecuencia del alelo E3: 0.8
Frecuencia del Alelo E4: 0.1
Se asume HWE
Riesgo Relativo de los genotipos E2E2 y E2E3: 0.8 (Se asume que ambos genotipos tienen el mismo RR)
Riesgo Relativo de los genotipos E3E4 y E2E4: 3 (Se asume que ambos genotipos tienen el mismo RR)
Los valores de Riesgo Relativo se basan en la tabla 1 del artículo “Systematic meta-analyses of Alzheimer disease genetic association studies: the AlzGene database”
ESTA SIMULACIÓN ES REPRODUCIBLE
Este script es completamente reproducible. Para obtener el cálculo de poder estadístico para otras condiciones, modificar los parámetros de la parte inferior y copiar todo el código en R.
Si el codigo no funciona, es posible que falten instalar algunos paquetes, de ser el caso, copiar lo siguiente en una ventana de R antes de empezar la simulación:
install.packages(“adegenet”) El paquete Rassoc ya no se encuentra disponible en el repositorio CRAN. Descargarlo manualmente e instalarlo.
# Parámetros de la simulación: Modificar según se desee
options(warn = F)
rr44_vector <- c(5, 7, 8, 9 , 10) ## Valores de Riesgo Relativo del genotipo E4E4 que serán ploteados
sample_size = c(50, 100, 150, 200, 250, 300) ## Tamaños muestrales que serán ploteados
preval = 0.0345 ## prevalencia
prop_casos = 0.5 ## proporción (número de casos)/(número total de individuos)
nreplicas = 1e4
## Número de réplicas empleadas en la simulación para calcular el poder estadístico
## Para 1000 replicas la desviación estándar máxima es 0.0158 (1.58%)
## Para 10000 replicas la desviación estándar máxima es 0.005 (0.5%)
p2 = 0.1 #Frecuencia del alelo E2
p3 = 0.8 #Frecuencia del alelo E3
p4 = 0.1 #Frecuencia del alelo E4
rr_22 = 0.8 ## Riesgo relativo de los genotipos E2E2 y E2E3 (Se asume que ambos genotipos tienen el mismo RR)
rr_23 = 0.8
rr_24 = 3
rr_34 = 3 ## Riesgo relativo de los genotipos E3E4 y E2E4 (Se asume que ambos genotipos tienen el mismo RR)
## Se establece el RR de E3E3 = 1 (es decir, como punto de referencia)INICIO DEL CÓDIGO QUE EJECUTA LA SIMULACIÓN
## NO MODIFICAR NADA DEBAJO DE ESTA LINEA
#########################################
require(adegenet)
require(Rassoc)
require(adegenet)
perform_assoc <- function(freq_e2 = 0.15,
freq_e3 = 0.8,
freq_e4 = 0.05,
rr22 = 0.7,
rr23 = 0.7,
rr24 = 2,
rr34 = 2,
rr44 = rr44_vector,
prevalencia = preval,
ncasos = 1e3,
ncontroles = 1e3,
met = metodo){
# freq_e2 = 0.1
# freq_e3 = 0.8
# freq_e4 = 0.1
#
# rr22 = 1
# rr23 = 1
# rr24 = 1
# rr34 = 1
# rr44 = c(2,4,5, 10)
# prevalencia = 0.0345
# ncasos = 1e2
# ncontroles = 1e2
# rr33 = 1
#
#
rr33 = 1
rr <- cbind(rr22, rr23, rr33, rr24, rr34, rr44) ## riesgo relativo
fgenotipicas <-c(freq_e2^2, 2*freq_e2*freq_e3,
freq_e3^2,
2*freq_e4* freq_e2, 2*freq_e4 * freq_e3,
freq_e4^2) ## frecuencias genotipicas
teta <- rr %*%fgenotipicas
w = prevalencia*rr/as.numeric(teta) ## Riesgo Absoluto
frecuencias_en_casos <- t(t(w)*fgenotipicas)/prevalencia ## frecuencias genotipicas relativas en los individuos afectados
#rowSums(frecuencias_en_casos) ## Sanity check: debe sumar 1
frecuencias_en_controles <- t(t((1-w))*fgenotipicas)/(1-prevalencia) ## frecuencias genotipicas relativas en los individuos controles
if(met == "logistic"){
# genotipos <- c("e22", "e23", "base_e33", "e24", "e34", "e44")
genotipos <- c(0, 0, 0, 1, 1, 2) ## "e22", "e23", "base_e33", "e24", "e34", "e44"
## Cada columna corresponde a un RR44 diferente:
casos <- apply(frecuencias_en_casos, MARGIN = 1,
function(x) sample(genotipos, size = ncasos, replace = T, prob =x))
controles <- apply(frecuencias_en_controles, MARGIN = 1,
function(x) sample(genotipos, size = ncontroles, replace = T, prob = x))
df <- data.frame(rbind(casos, controles))
status = c(rep(1, ncasos), rep(0, ncontroles))
apply(df, 2, function(x)summary(glm(status ~ x, family = "binomial"))$coefficients[2,4])
} else if (met == "allelic") {
genotipos <- c(0, 0, 0, 1, 1, 2)
casos <- apply(frecuencias_en_casos, MARGIN = 1,
function(x) sample(genotipos, size = ncasos, replace = T, prob =x))
controles <- apply(frecuencias_en_controles, MARGIN = 1,
function(x) sample(genotipos, size = ncontroles, replace = T, prob = x))
df <- data.frame(rbind(casos, controles))
status = c(rep(1, ncasos), rep(0, ncontroles))
apply(df, 2, function(x) chisq.test(table(factor(x, levels = c(0,1,2)),
status)[1:2,] + rbind(c(0,0),table(factor(x, levels = c(0,1,2)), status)[3,]*2))$p.val)
} else {
genotipos <- c("e22", "e23", "e33", "e24", "e34", "e44")
casos <- apply(frecuencias_en_casos, MARGIN = 1,
function(x) sample(genotipos, size = ncasos, replace = T, prob = x))
controles <- apply(frecuencias_en_controles, MARGIN = 1,
function(x) sample(genotipos, size = ncontroles, replace = T, prob = x))
grupos <- vector("list", length(rr44))
for(k in 1:length(rr44)){
afectados <- as.numeric(table(factor(casos[,k], levels = genotipos)))
sanos <- as.numeric(table(factor(controles[,k], levels = genotipos)))
grupos[[k]] <- data.frame(afectados, sanos)
}
if(met == "chi2"){
sapply(grupos, function(x) chisq.test(rbind(colSums(x[1:3,]), colSums(x[4:6,])), correct = F)$p.value)
} else if(met=="catt"){
sapply(grupos, function(x) CATT(cbind(rowSums(t(x)[,1:3]),rowSums(t(x)[,4:5]), t(x)[,6]))$p.value)
} else {
print ("Método de asociación incorrecto. Elegir entre 'chi2' para la prueba de Chi-cuadrado, 'catt' para la prueba de Cochran-Armitage Trend test', y 'logistic' para regresión logística")
}
}
}RESULTADOS: CHI-CUADRADO
options(warn = F)
metodo = "chi2"
power_matrix <- data.frame(sapply(sample_size, function(x) (rowMeans(replicate(nreplicas,
perform_assoc(freq_e2 = p2,
freq_e3 = p3,
freq_e4 = p4,
rr22 = rr_22,
rr23 = rr_23,
rr24 = rr_24,
rr34 = rr_34,
rr44 = rr44_vector,
met = metodo,
prevalencia = preval,
ncasos = round(x * prop_casos ),
ncontroles = round(x *(1-prop_casos) ) )) < 0.05))))
colnames(power_matrix) <- paste0("n=", sample_size)
rownames(power_matrix) <- paste0("RR44=", rr44_vector)
colores <- adjustcolor(funky(8)[1:length(rr44_vector)], alpha.f = 0.7)
plot(sample_size, power_matrix[1,], type = "b", col = colores[1], lwd = 3,
main = "Poder estadístico de la prueba Chi-cuadrado",
xlab = "Tamaño muestral", ylab = "Poder estadístico", cex.lab = 1.2, ylim = c(0,1))
for(i in 2:length(rr44_vector)){
points(sample_size, power_matrix[i,], type = "b", col = colores[i], lwd = 3)
}
grid()
leg_texto = paste0("RR E4E4 = ", rr44_vector)
legend("bottomright", legend = leg_texto, col = colores, pch = 16, bty = "n", pt.cex = 1.5)
| n=50 | n=100 | n=150 | n=200 | n=250 | n=300 | |
|---|---|---|---|---|---|---|
| RR44=5 | 0.5077 | 0.8016 | 0.9246 | 0.9763 | 0.9919 | 0.9979 |
| RR44=7 | 0.5252 | 0.8162 | 0.9392 | 0.9829 | 0.9961 | 0.9993 |
| RR44=8 | 0.5403 | 0.8188 | 0.9453 | 0.9849 | 0.9959 | 0.9984 |
| RR44=9 | 0.5521 | 0.8489 | 0.9516 | 0.9869 | 0.9967 | 0.9994 |
| RR44=10 | 0.5637 | 0.8508 | 0.9541 | 0.9891 | 0.9965 | 0.9994 |
RESULTADOS: COCHRAN-ARMITAGE TREND TEST
options(warn = F)
metodo = "catt"
power_matrix <- data.frame(sapply(sample_size, function(x) (rowMeans(replicate(nreplicas,
perform_assoc(freq_e2 = p2,
freq_e3 = p3,
freq_e4 = p4,
rr22 = rr_22,
rr23 = rr_23,
rr24 = rr_24,
rr34 = rr_34,
rr44 = rr44_vector,
met = metodo,
prevalencia = preval,
ncasos = round(x * prop_casos ),
ncontroles = round(x *(1-prop_casos) ) )) < 0.05))))
colnames(power_matrix) <- paste0("n=", sample_size)
rownames(power_matrix) <- paste0("RR44=", rr44_vector)
colores <- adjustcolor(funky(8)[1:length(rr44_vector)], alpha.f = 0.7)
plot(sample_size, power_matrix[1,], type = "b", col = colores[1], lwd = 3,
main = "Poder estadístico de la prueba Cochran-Armitage Trend Test",
xlab = "Tamaño muestral", ylab = "Poder estadístico", cex.lab = 1.2, ylim = c(0,1))
for(i in 2:length(rr44_vector)){
points(sample_size, power_matrix[i,], type = "b", col = colores[i], lwd = 3)
}
grid()
leg_texto = paste0("RR E4E4 = ", rr44_vector)
legend("bottomright", legend = leg_texto, col = colores, pch = 16, bty = "n", pt.cex = 1.5)
| n=50 | n=100 | n=150 | n=200 | n=250 | n=300 | |
|---|---|---|---|---|---|---|
| RR44=5 | 0.4039 | 0.7686 | 0.9160 | 0.9738 | 0.9920 | 0.9972 |
| RR44=7 | 0.4476 | 0.8171 | 0.9410 | 0.9853 | 0.9957 | 0.9994 |
| RR44=8 | 0.4663 | 0.8342 | 0.9546 | 0.9900 | 0.9964 | 0.9996 |
| RR44=9 | 0.4853 | 0.8593 | 0.9643 | 0.9936 | 0.9982 | 0.9997 |
| RR44=10 | 0.5263 | 0.8663 | 0.9690 | 0.9938 | 0.9988 | 0.9997 |
RESULTADOS: REGRESIÓN LOGÍSTICA
options(warn = F)
metodo = "logistic"
power_matrix <- data.frame(sapply(sample_size, function(x) (rowMeans(replicate(nreplicas,
perform_assoc(freq_e2 = p2,
freq_e3 = p3,
freq_e4 = p4,
rr22 = rr_22,
rr23 = rr_23,
rr24 = rr_24,
rr34 = rr_34,
rr44 = rr44_vector,
met = metodo,
prevalencia = preval,
ncasos = round(x * prop_casos ),
ncontroles = round(x *(1-prop_casos) ) )) < 0.05))))
colnames(power_matrix) <- paste0("n=", sample_size)
rownames(power_matrix) <- paste0("RR44=", rr44_vector)
colores <- adjustcolor(funky(8)[1:length(rr44_vector)], alpha.f = 0.7)
plot(sample_size, power_matrix[1,], type = "b", col = colores[1], lwd = 3,
main = "Poder estadístico de la regresión logística",
xlab = "Tamaño muestral", ylab = "Poder estadístico", cex.lab = 1.2, ylim = c(0,1))
for(i in 2:length(rr44_vector)){
points(sample_size, power_matrix[i,], type = "b", col = colores[i], lwd = 3)
}
grid()
leg_texto = paste0("RR E4E4 = ", rr44_vector)
legend("bottomright", legend = leg_texto, col = colores, pch = 16, bty = "n", pt.cex = 1.5)
| n=50 | n=100 | n=150 | n=200 | n=250 | n=300 | |
|---|---|---|---|---|---|---|
| RR44=5 | 0.4449 | 0.7849 | 0.9240 | 0.9744 | 0.9929 | 0.9970 |
| RR44=7 | 0.4744 | 0.8239 | 0.9446 | 0.9847 | 0.9969 | 0.9992 |
| RR44=8 | 0.5044 | 0.8420 | 0.9538 | 0.9896 | 0.9967 | 0.9992 |
| RR44=9 | 0.5126 | 0.8621 | 0.9637 | 0.9932 | 0.9982 | 0.9997 |
| RR44=10 | 0.5385 | 0.8784 | 0.9673 | 0.9940 | 0.9991 | 0.9998 |
RESULTADOS: ASOCIACIÓN ALÉLICA
options(warn = F)
metodo = "allelic"
power_matrix <- data.frame(sapply(sample_size, function(x) (rowMeans(replicate(nreplicas,
perform_assoc(freq_e2 = p2,
freq_e3 = p3,
freq_e4 = p4,
rr22 = rr_22,
rr23 = rr_23,
rr24 = rr_24,
rr34 = rr_34,
rr44 = rr44_vector,
met = metodo,
prevalencia = preval,
ncasos = round(x * prop_casos ),
ncontroles = round(x *(1-prop_casos) ) )) < 0.05))))
colnames(power_matrix) <- paste0("n=", sample_size)
rownames(power_matrix) <- paste0("RR44=", rr44_vector)
colores <- adjustcolor(funky(8)[1:length(rr44_vector)], alpha.f = 0.7)
plot(sample_size, power_matrix[1,], type = "b", col = colores[1], lwd = 3,
main = "Poder estadístico de la asociación alélica",
xlab = "Tamaño muestral", ylab = "Poder estadístico", cex.lab = 1.2, ylim = c(0,1))
for(i in 2:length(rr44_vector)){
points(sample_size, power_matrix[i,], type = "b", col = colores[i], lwd = 3)
}
grid()
leg_texto = paste0("RR E4E4 = ", rr44_vector)
legend("bottomright", legend = leg_texto, col = colores, pch = 16, bty = "n", pt.cex = 1.5)
| n=50 | n=100 | n=150 | n=200 | n=250 | n=300 | |
|---|---|---|---|---|---|---|
| RR44=5 | 0.4081 | 0.7628 | 0.9171 | 0.9752 | 0.9915 | 0.9974 |
| RR44=7 | 0.4555 | 0.8032 | 0.9420 | 0.9833 | 0.9956 | 0.9994 |
| RR44=8 | 0.4831 | 0.8214 | 0.9513 | 0.9883 | 0.9979 | 0.9993 |
| RR44=9 | 0.5073 | 0.8368 | 0.9617 | 0.9895 | 0.9982 | 0.9998 |
| RR44=10 | 0.5263 | 0.8656 | 0.9678 | 0.9922 | 0.9989 | 1.0000 |
CONFIRMACIÓN DE QUE LA SIMULACIÓN FUNCIONA CORRECTAMENTE
Notar que para ninguno de los métodos de asociación existe mucha diferencia en el poder estadístico para diferentes valores de RR del genotipo E4E4. Ello se debe a que el poder estadístico está gobernado principalmente por el genotipo más frecuente; es decir, el genotipo E3E4.
Para comprobar que el enunciado anterior es correcto, basta con repetir la simulación pero con valores de RR igual a 1 para todos los genotipos excepto el E4E4. Además, se incluirá como referencia el model nulo, en el que todos los genotipos (incluyendo E3E4 y E4E4) tienen un RR = 1. El poder estadístico aparente en el último escenario debe ser de 0.05, que corresponde al alfa seleccionado. Solo se muestra la comprobación para el método Chi-cuadrado.
rr44_vector <- c(1, 5, 7, 8, 9 , 10)
sample_size <- c(100, 300, 500, 800, 1000, 1300)
metodo = "chi2"
power_matrix <- data.frame(sapply(sample_size, function(x) (rowMeans(replicate(nreplicas,
perform_assoc(freq_e2 = p2,
freq_e3 = p3,
freq_e4 = p4,
rr22 = 1,
rr23 = 1,
rr24 = 1,
rr34 = 1,
rr44 = rr44_vector,
met = metodo,
prevalencia = preval,
ncasos = round(x * prop_casos ),
ncontroles = round(x *(1-prop_casos) ) )) < 0.05))))
colnames(power_matrix) <- paste0("n=", sample_size)
rownames(power_matrix) <- paste0("RR44=", rr44_vector)
colores <- adjustcolor(funky(8)[1:length(rr44_vector)], alpha.f = 0.7)
plot(sample_size, power_matrix[1,], type = "b", col = colores[1], lwd = 3,
xlab = "Tamaño muestral", ylab = "Poder estadístico", cex.lab = 1.2, ylim = c(0,1))
for(i in 2:length(rr44_vector)){
points(sample_size, power_matrix[i,], type = "b", col = colores[i], lwd = 3)
}
grid()
leg_texto = c(paste0("RR E4E4 = ", rr44_vector), "alpha = 0.05")
legend("topleft", legend = leg_texto, col = c(colores, adjustcolor("gray60", alpha.f = 0.5)),
pch = 16, bty = "n", pt.cex = 1.5)
abline(h = 0.05, lwd = 2, lty = 2, col = adjustcolor("gray60", alpha.f = 0.5))
kable(power_matrix, caption = "Poder estadístico de la prueba Chi-cuadrado")| n=100 | n=300 | n=500 | n=800 | n=1000 | n=1300 | |
|---|---|---|---|---|---|---|
| RR44=1 | 0.0529 | 0.0520 | 0.0525 | 0.0515 | 0.0566 | 0.0509 |
| RR44=5 | 0.0694 | 0.1098 | 0.1422 | 0.1988 | 0.2425 | 0.3002 |
| RR44=7 | 0.0871 | 0.1689 | 0.2640 | 0.3685 | 0.4463 | 0.5416 |
| RR44=8 | 0.0984 | 0.2203 | 0.3216 | 0.4715 | 0.5609 | 0.6745 |
| RR44=9 | 0.1180 | 0.2583 | 0.3894 | 0.5597 | 0.6576 | 0.7717 |
| RR44=10 | 0.1350 | 0.3054 | 0.4612 | 0.6557 | 0.7519 | 0.8535 |
Se demuestra que la frecuencia de errores tipo I en el modelo nulo (el RR de todos los genotipos es igual a 1) es igual al alfa nominal de 0.05; Por lo tanto, la simulación es correcta y se confirma que en las 4 gráficas iniciales el poder estadístico es semejante en todos los escenarios debido a que este está gobernado por el RR del genotipo E3E4 por ser este mucho más frecuente.