famuvie
commented
2 years ago Dear Gregor,
Thanks for your interesting question, and for your clearly stated reproducible example!
You compare results of accuracy of predicted genetic values from breedR and rrBLUP out of simulated data and show that for about 23 % of the cases, the results from breedR collapse down to about 0, whereas the results from rrBLUP remain around the expected value.
I think the issue boils down to numerical instability, where rrBLUP stands more robustly thanks to (I'm guessing) its built-in shrinkage capabilities that prevent variance estimates to collapse.
However, you are measuring accuracy via the correlation between predicted and true genetic values, which is not necessarily a good measure of predictive performance. As I show below, breedR and rrBLUP provide very similar predictions, even in the cases where accuracies are very different.
To demonstrate this, let me examine more in detail the results of one of those badly behaving cases:
# Packages ----------------------------------------------------------------
library(breedR)
#> Loading required package: sp
library(rrBLUP)
# Functions ---------------------------------------------------------------
draw_M <- function(rows, cols) {
values <- 1 - 2 * (runif(rows * cols) < 0.5)
return(matrix(values, rows, cols))
}
simulate_data <- function(n_obs, n_marks, h2) {
M <- draw_M(n_obs, n_marks)
u <- rnorm(n_marks)
g <- as.vector(M %*% u)
y <- g + rnorm(n_obs, mean = 0, sd = sqrt((1 - h2) / h2 * var(g)))
structure(
data.frame(g, y),
G = A.mat(M)
)
}
# A badly behaved case ----------------------------------------------------
set.seed(20220403)
baddat <- simulate_data(n_obs = 200, n_marks = 1e3, h2 = 0.5)
ans <- mixed.solve(baddat$y, K = attr(baddat, 'G'))
(accuracy <- cor(baddat$g, ans$u))
#> [1] 0.7073303
#predict breeding values with breedR
model <- remlf90(
fixed = y ~ 1,
generic = list(gen = list(diag(nrow(baddat)), attr(baddat, 'G'))),
method = 'ai',
data = baddat
)
#> Using default initial variances given by default_initial_variance()
#> See ?breedR.getOption.
(accuracy_b <- cor(baddat$g, model$ranef$gen[[1]]$value))
#> [1] -0.04447799
preds <- baddat |>
transform(
pred_rr = ans$u,
pred_breedR = model$ranef$gen[[1]]$value
)
pairs(data.matrix(preds), xlim = c(-150, 150), ylim = c(-150, 150))
## The variance of the genetic effect was estimated at essentially 0,
## collapsing all estimated genetic values at 0 (except for 1 that deviates a
## little bit, for some numerical artefact or idiosyncrasy).
summary(model)
#> Formula: y ~ 0 + Intercept
#> Data: baddat
#> AIC BIC logLik
#> 2038 2045 -1017
#>
#> Parameters of special components:
#>
#>
#> Variance components:
#> Estimated variances S.E.
#> gen 1.299e-04 1.306e-05
#> Residual 1.720e+03 1.720e+02
#>
#> Fixed effects:
#> value s.e.
#> Intercept -3.4838 2.94Created on 2022-04-04 by the reprex package (v2.0.1)
As you see, the accuracies from breedR and rrBLUP are dramatically different (near 0 and 70 % respectively). Yet, the predicted genetic values are essentially the same (see bottom-right panels of the pairs-plot, where predictions from both programs are basically collapsed over the same spot, when compared to the variation of the true genetic values).
breedR's estimate of the variance has collapsed to zero, and all the predicted values are essentially 0, except for one slightly different value, which honestly surprised me, but must have something to do with the idiosyncrasies of the underlying algorithm. This is responsible of the violent bad accuracy (= correlation). But, as you can see, it is simply an artefact of the very bad-conditioned data, which rrBLUP handles certainly more elegantly.
GregorDall
Dear Facundo,
we are trying to implement GBLUP / rrBLUP using breedR but seem to have stumbled upon an issue. In some cases breedR produces unreasonable results and also takes comparably longer to produce an output. Here is an example comparing the results obained by rrBLUP and breedR. Is there any way you can explain this behaviour or to avoid it? I have also tried using "em", but the results are the same
Cheers Gregor
` library(breedR) library(rrBLUP)
result_all <- NULL
for(i in 1:100){
cat(paste("Run",i,"at",date(),"\n"))
random population of 200 lines with 1000 markers
M <- matrix(rep(0,200*1000),200,1000) for (i in 1:200) { M[i,] <- ifelse(runif(1000)<0.5,-1,1) }
random phenotypes
u <- rnorm(1000) g <- as.vector(crossprod(t(M),u)) h2 <- 0.5 #heritability y <- g + rnorm(200,mean=0,sd=sqrt((1-h2)/h2*var(g)))
predict breeding values with rrBLUP
G <- A.mat(M) ans <- mixed.solve(y,K=G) accuracy <- cor(g,ans$u)
predict breeding values with breedR
model <- remlf90(fixed = y ~ 1, generic = list(gen = list(diag(200),G)), method = 'ai', data = as.data.frame(y)) accuracy_b <- cor(g,model$ranef$gen[[1]]$value)
result <- data.frame(accuracy,accuracy_b) result_all <- rbind(result_all,result)
} summary(result_all)#expectation r~sqrt(0.5)~0.7071068 length(which(result_all$accuracy>result_all$accuracy_b)) length(which(result_all$accuracy<0.5)) length(which(result_all$accuracy_b<0.5))
par(mfrow=c(1,3)) hist(result_all$accuracy,main="rrBLUP",xlim=c(-0.25,1),xlab="prediction accuracy") hist(result_all$accuracy_b,main="breedR",xlim=c(-0.25,1),xlab="prediction accuracy") plot(result_all$accuracy,result_all$accuracy_b,xlim=c(-0.25,1),ylim=c(-0.25,1), xlab="rrBLUP",ylab="breedR") abline(a=0,b=1,col="red") `