mc.se: Count error for h2 and corr

In yzhlinscau/AAfun: ASReml-R Added functions

Description

This function counts standard error(se) for heritability(h2) and corr value for MCMCglmm package.

Usage

mc.se(object = NULL, Nmc = NULL, confinterval = NULL, lv = NULL, 
        uv = NULL, n = NULL, conf.level = NULL, sigf = NULL)

Arguments

object	MCMCglmm model results or h2/corr formula results
Nmc	Use MCMCglmm result directly (T) or not (F)
confinterval	Confidence interval for heritability or corr
lv	Lower value from confidence interval
uv	Upper value from confidence interval
n	Total number of aim trait observation
conf.level	Confidence level
sigf	Output significent level(T) or not(F,default)

Value

Nmc	Default is T to use MCMCglmm results directly, F for not.
n	Total number of aim trait observation,default value is 1000.
conf.level	Confidence level, default value is 0.95.
sigf	Output significant level (T) or not (default F).

Author(s)

Yuanzhen Lin <[e-mail]>

References

Yuanzhen Lin. R & ASReml-R Statistics. China Forestry Publishing House. 2016

See Also

Website for instant update: yzhlin-asreml.ys168.com

Examples

library(AAfun)
library(MCMCglmm)
        
data(PrSpa)

df<-subset(PrSpa,Spacing=='3')

######## single trait model ######## 
m1.glmm <- MCMCglmm(h5 ~ 1 + Rep, 
                    random = ~ Fam, pr = TRUE,
                    family = 'gaussian',
                    data = df)

h2.glmm<-4*m1.glmm$VCV[,'Fam']/(m1.glmm$VCV[,'Fam']+m1.glmm$VCV[,'units'])
posterior.mode(h2.glmm)

mc.se(h2.glmm)

confinterval<-HPDinterval(h2.glmm)
mc.se(confinterval=confinterval,Nmc=F)
mc.se(confinterval=confinterval,Nmc=F,n=559,conf.level=0.95)

### Second method has the same result ###
lv<-HPDinterval(h2.glmm)[1]
uv<-HPDinterval(h2.glmm)[2]
mc.se(lv=lv,uv=uv,Nmc=F)

######## bi-trait model ########

df$dj<-1000*df$dj
phen.var<-matrix(c(var(df$dj,na.rm=TRUE),0,0,
                   var(df$h5,na.rm=TRUE)),2,2)

prior<-list(G=list(G1=list(V=phen.var/2,n=2)),
            R=list(V=phen.var/2,n=2))

set.seed(1234)
m2.glmm <- MCMCglmm(cbind(dj,h5)~ trait-1+trait:Rep, random=~ us(trait):Fam, 
                    rcov=~ us(trait):units, data=df, family=c("gaussian", "gaussian"), 
                    nitt=130000,thin=100,burnin=30000,
                    prior=prior,verbose=FALSE,pr=TRUE)

posterior.mode(m2.glmm$VCV)
HPDinterval(m2.glmm$VCV)
#### count se for variance component
mc.se(m2.glmm$VCV)

#### count se for fixed and randomed effects
# mc.se(m2.glmm$Sol)
posterior.mode(m2.glmm$Sol)[c(1:5,40:45,80:85)]
mc.se(m2.glmm$Sol)[c(1:5,40:45,80:85),]

#### count se for heritability
A.h2.glmm<-4*m2.glmm$VCV[,'dj:dj.Fam']/(m2.glmm$VCV[,'dj:dj.Fam']+m2.glmm$VCV[,'dj:dj.units'])
posterior.mode(A.h2.glmm)

mc.se(A.h2.glmm)

confinterval<-HPDinterval(A.h2.glmm)
mc.se(confinterval=confinterval,Nmc=F)
mc.se(confinterval=confinterval,n=559,conf.level=0.95,Nmc=F)

### Second method has the same result for h2 ###
lv<-HPDinterval(A.h2.glmm)[,1]
uv<-HPDinterval(A.h2.glmm)[,2]
mc.se(lv=lv,uv=uv,Nmc=F)

#### count se for corr
gCorr.glmm<-m2.glmm$VCV[,'h5:dj.Fam']/sqrt(m2.glmm$VCV[,'dj:dj.Fam']*m2.glmm$VCV[,'h5:h5.Fam'])
mc.se(gCorr.glmm,sigf=T)

yzhlinscau/AAfun documentation built on May 21, 2020, 2:19 p.m.