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 [email protected] 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
#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)