data(ratones)
# Estimating a W matrix, controlling for line and sex
model_formula = paste0("cbind(",
paste(names(ratones)[13:47], collapse = ", "),
") ~ SEX + LIN")
ratones_W_model = lm(model_formula, data = ratones)
W_matrix = CalculateMatrix(ratones_W_model)
# Estimating the divergence between the two direction of selection
delta_Z = colMeans(ratones[ratones$selection == "upwards", 13:47]) -
colMeans(ratones[ratones$selection == "downwards", 13:47])
# Reconstructing selection gradients with and without noise control
Beta = solve(W_matrix, delta_Z)
Beta_non_noise = solve(ExtendMatrix(W_matrix, ret.dim = 10)$ExtMat, delta_Z)
# Comparing the selection gradients to the observed divergence
Beta %*% delta_Z /(Norm(Beta) * Norm(delta_Z))
Beta_non_noise %*% delta_Z /(Norm(Beta_non_noise) * Norm(delta_Z))
Run the code above in your browser using DataLab