Poder estadistico - Rasgo cunatitativo

CALCULO DE LA POTENCIA ESTADISTICA MEDIANTE MÉTODO ANALITICO

Solo se incluyen los efectos aditivos en el modelo

N = 90
alpha = 0.05
H2 = seq(0.02, 0.16, 0.02)

threshold = qchisq(alpha, df = 1, lower.tail = FALSE)
poder = pchisq(threshold, df = 1, lower.tail = FALSE, ncp = N * H2)

par(new=T)
plot(H2, poder, type = "o", col = "gray60", lwd = 1, pch = 19, cex = 1,
     xlab = "Varianza genetica aditiva/Varianza total", ylab= "Potencia estadistica",
     main = "Potencia estadistica vs varianza genética aditiva debida al QTL", cex.lab = 1.2)
grid(lwd = 2)

CALCULO DE LA POTENCIA ESTADISTICA MEDIANTE SIMULACION

simulation <- function(nind, p, a, d, total.sd, p.masc, sex.effect, modelo = "aditivo"){
    q = 1-p
    alle = c(1, 0)
    allefreq = c(p, q)
    sexos = c(0,1)
    p.fem = 1-p.masc
    prob.sexo = c(p.fem, p.masc)
    data <- matrix(0, nind, 6)
    
    alfa = a + d*(q-p)
    varA = p * q * alfa^2
    varD = (2*p*q*d)^2
    varSex = p.masc * p.fem * sex.effect^2
    s = sqrt(total.sd^2 -( varA + varD + varSex))
    
    data <- matrix(sample(alle, size = 2*nind, replace = T, prob= allefreq), ncol = 2)
    score_ad <- rowSums(data)-1
    score_dom <- as.integer(rowSums(data) == 1)
    sexo <- sample(sexos, size = nind, replace = T, prob = prob.sexo)
    fenotipo <- score_ad*a + rnorm(nind, 0 , s) + sexo * sex.effect
    
    if(modelo =="aditivo"){
        pval <- summary(lm(fenotipo ~ score_ad))$coefficients[2,4]
    } else if( modelo == "dominante"){
        pval <- summary(lm(fenotipo ~ score_ad + score_dom))$coefficients[2,4]
    }
        
        pval < 0.05
}

## Check 1: funcion "simulation" funciona adecuadamente
# simulation(nind = 152, a = 1, d = 1 , p = 0.5, total.sd = 3.04, p.masc = 0.0974, sex.effect = 0.8)

mean.power <- function(nind, a, d , p, total.sd, p.masc, n.reps, sex.effect, mod = "aditivo") {
    reps <- replicate(n.reps, simulation(nind = nind, a = a, d = d, p = p,
                                         p.masc = p.masc , total.sd = total.sd,
                                         sex.effect = sex.effect, modelo = mod))
    return(mean(reps))
}

Parametros de la simulacion

a_vector = seq(0.3, 1.5, 0.3)
p_vector = seq(0.1, 0.9, 0.2)
total.sd = 3.03
n.reps = 1e4
nind = 145
sex.effect = 0.8
p.masc = 0.1 ##0.0974

Resultados: Rasgo con dominancia parcial, modelo ajustado = aditivo

power.matrix <- matrix(0 , nrow = length(a_vector), ncol = length(p_vector))

for( i in 1:length(p_vector)){
    power.matrix[,i] <- sapply(a_vector, function(a_vector){ mean.power(nind = nind, a = a_vector, d = a_vector/3, 
                                                                        total.sd = total.sd,  p = p_vector[i],
                                                                        n.reps = n.reps, p.masc = p.masc,
                                                                        sex.effect = sex.effect, mod = "aditivo")})
}

## Cada columna de power.matrix contiene el poder estadistico para un valor de p,
# que coincide con cada columna de la matriz h2.matrix a continuacion:

## Crear una matrix tal que:
# Numero de filas = Numero de valores de a (h2 )
# Numero de columnas = Numero de valores de p
# Por cada valor de p se hara una curva diferente
# por cada valor de "r2" (cada valor de a) corresponde un punto en el eje x

h2.matrix = matrix(0, nrow = length(a_vector), ncol = length(p_vector))
p.fem = 1-p.masc
for( i in 1:length(p_vector)){
    q = 1 - p_vector[i]
    d = a_vector/3
    
    vara <- 2 * p_vector[i] * q *((a_vector + (q-p_vector[i])*d)^2)
    vart <- total.sd^2
    h2.matrix[,i] <- round(vara/vart, 4)
}

## A continuacion se graficara una curva de poder estadistico para cada p: Notar que
# el poder estadistico no solo depende de la porcion de la varianza explicada por el modelo
# sino tambien por la frecuencia alelica

len = length(p_vector)-1
cols = adjustcolor(c("orangered4", "lightblue","dodgerblue4", "ivory4", "darkgreen"), alpha.f = 0.50)

leyenda = vector("character", length(p_vector))
for(i in 1:length(p_vector)){
    leyenda[i] = paste("p = ", p_vector[i], sep = "")
}


plot(h2.matrix[,length(p_vector)], power.matrix[,length(p_vector)],
     type = "b", pch = 20, col = cols[length(p_vector)], lwd = 2,  
     ylim = c(0, 1), xlim = c(0, max(h2.matrix)*1.1),
     xlab = "Varianza genética aditiva/Varianza total", 
     ylab = "Potencia estadísstica",
     main = "Potencia estadística vs Varianza genética aditiva") 
grid(lwd = 2)
for(i in 1:len){
    lines(h2.matrix[,i], power.matrix[,i], type = "b", lwd = 2, pch = 20, col = cols[i])
}
legend("bottomright", legend = leyenda, bty = "n", lty = 7, lwd=c(3,3,3,3,3),col=cols)


par(new=T)
plot(H2, poder, type = "o", col = "gray70", lwd = 1, pch = 19, cex = 1,
     cex.lab = 1.2, xlab = "", ylab = "", axes = F, lty = 2)

Resultados: Rasgo con dominancia parcial, modelo ajustado = efectos aditivo y de dominancia

Resultados: Rasgo con dominancia completa, modelo ajustado = aditivo

Resultados: Rasgo con dominancia completa, modelo ajustado = efectos aditivo y de dominancia