Reliability estimate

[makgeha]

Gustavo, Is there any built in function to calculate a reliability estimate for each of the predicted values in the BGLR code? (this is similar to computing accuracy of prediction in Pedigree BLUP sqrt(1-(std_err_pred_i)^2/(diag(Aii)*s2a)) Thanks, Mak

[gdlc]

fm$SD.yHat

will give you the (estimated) posterior standard deviation...which is the bayesian cournterpart of the PEV

On Tue, Aug 9, 2016 at 2:47 PM, makgeha notifications@github.com wrote:

Gustavo, Is there any built in function to calculate a reliability estimate for each of the predicted values in the BGLR code? (this is similar to computing accuracy of prediction in Pedigree BLUP sqrt(1-(std_err_pred_i)^2/(diag(Aii)*s2a)) Thanks, Mak

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/gdlc/BGLR-R/issues/9, or mute the thread https://github.com/notifications/unsubscribe-auth/AE0NonlvBCbBWB4BkJ37AzMkgy5o8t5iks5qeMs_gaJpZM4JgZYJ .

[makgeha]

Gustavo, Thanks for your prompt reply. So for each model used (BRR, BL, BayesB etc...) I can get the SD.yhat for each individual. Where would be the estimated value of the additive genetic variance to plug into the accuracy function? Thanks, Mak

[cheuerde]

I am replying here some years later just because somebody fell over this issue as well and I would like to address the question regarding the accuracy (reliability) of breeding values.

library(BGLR)

# example dataset
data(wheat)

# pedigree A
A <- wheat.A

# response
y <- wheat.Y[,1]

# genotypes
M <- wheat.X

# G matrix
vars <- apply(M, 2, var)
G <- tcrossprod(scale(M, scale = FALSE)) / sum(vars)

# add something to diagional to stabilize things
diag(G) <- diag(G) + 0.01

# using L as design matrix will lead to an equivalent animal model, which is just
# nicer to handle in BGLR
L <- t(chol(G))

# fixed effects (just intercept)
X <- array(1, dim = c(length(y), 1))

# run animal model
niter = 20000
burnin = 5000
mod <- BGLR(
        y = y,
        response_type = "gaussian",
        ETA = list(
               list(X = X, model = "FIXED"),
               # list(K = G, model = "RKHS", saveEffects = TRUE)
               list(X = L, model = "BRR", saveEffects = TRUE)
               ),
        nIter = niter,
        burnIn = burnin,
        thin = 1,
        saveAt = "myModel",
        verbose = FALSE
)

# get the posterior (only stored in file)
uFile <- mod$ETA[[2]]$NamefileOut
uDat <- scan(uFile, what = numeric(), sep = "\n")

eDat <- scan("myModelvarE.dat",  what = numeric(), sep = "\n")

postit <- (burnin + 1) : niter

varU <- uDat[postit]
varE <- eDat[postit]
h2 <- varU / (varU + varE)

# the heritability and its standard deviation (~ standard error)
h2Estimate = mean(h2)
h2SD = sd(h2)

# now reliabilities. 
# we specified to store the posterior samples of the random effects with the option "saveEffects = TRUE".
# read those samples in now
u = readBinMat('myModelETA_2_b.bin')

# the breeding values (posterior)
g = L %*% t(u)

# PEV
PEV = apply(g, 1, var)

# reliabilities
vU <- mean(varU)
REL = 1 - (PEV / vU)