Poder estadistico - Rasgo cunatitativo

Statistical power through analytical method:

Solo se incluyen los efectos aditivos en el modelo

N = 140
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=F)
plot(H2, poder, type = "o", col = "gray60", lwd = 1, pch = 19, cex = 1,
     xlab = "Additive genetic variance/total variance", ylab= " Power",
     main = "Statistical power", cex.lab = 1.2)
grid()

Statistical power through simulations:

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
}

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))
}

Simulation parameters

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

Results: Locus with partial dominance, Adjusted model = additive

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")})
}

## Each column of power.matrix contains the power for a certain allele frequency "p",
# corresponding with each column of the matrix h2.matrix up next:

## Create a matrix where: 
# Number of rows = number of h2 values
# Number of columns = number of values of "p" (allele frequency)
# For each allele frequency, make a curve 

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)
}
## Plotting:

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 = "Additive genetic variance/Total variance", 
     ylab = "Power",
     main = "Statistical power vs additive genetic variance") 
grid()
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)

Results: Locus with partial dominance, adjusted model = additive + dominance effects

Results: locus with complete dominance, adjusted model = additive

Results: locus with complete dominance, adjusted model = additive and dominance effects