DIC 기반 이형질자체 유전 모델: 가산형 vs 우성형

두-알레일 모델과 DIC 비교 두-알레일 모델은 알레일 A1과 A2, 그리고 형질형식 A1A1, A1A2, A2A2 를 갖습니다. 이 모델은 additive effect만을 사용해(제한 모델) 또는 additive + dominance effect를 사용해(전체 모델) 적합할 수 있습니다. 전체 모델은 제한 모델보다 파라미터가 하나 더 많습니다. 최근에 이 모델을 사용해 DIC 결과를 얻었는데, 전체 모델이 파라미터가 적은 것처럼 보이는 이상한 결과가 나왔습니다. 아래는 R과 BUGS를 사용해 시뮬레이션한 예제로, 기대되는 DIC 결과를 보여줍니다. 필요한 패키지와 함수 library(package="e1071") library(package="doBy") library(package="R2WinBUGS") rvalue { if(!is.matrix(value)) value ## Frequencies k q P Q H ## Allocate matrix of arbitrary values by loci x y ## Simulate arbitrary values for(i in 1:k) { x[, i] y[, i] } x ## Sum all values y ## Add mean and normal deviate (error) y return(as.data.frame(cbind(y, x))) } prepairData { ret ret$n ret$y if(!dominance) { ret$x } else { ret$x } ret } 시뮬레이션 데이터 생성 ### SIMULATE DATA###----------------------------------------------------------------------- ## --- Purely additive --- set.seed(seed=19791123) ## Model parameters n mu sigma value p ## SimulatepodatkiA par(mfrow=c(2, 1)) hist(podatkiA$y, xlab="Phenotype") plot(podatkiA$y ~ jitter(podatkiA$V2), xlab="Genotype", ylab="Phenotype") tmpA1 abline(coef(tmpA1), lwd="2", col="blue") tmpA2 points(tmpA2$y.mean ~ tmpA2$V2, pch=19, col="red") ## --- Additive + Dominance --- set.seed(seed=19791123) ## Model parameters n mu sigma value p ## SimulatepodatkiD par(mfrow=c(2, 1)) hist(podatkiD$y, xlab="Phenotype") plot(podatkiD$y ~ jitter(podatkiD$V2), xlab="Genotype", ylab="Phenotype") tmpD1 abline(coef(tmpD1), lwd="2", col="blue") tmpD2 points(tmpD2$y.mean ~ tmpD2$V2, pch=19, col="red") DIC을 이용한 모델 파라미터와 추정 ### ESTIMATE MODEL PARAMETERS AND DIC###----------------------------------------------------------------------- modelA { ## --- Additive model (regression on gene content) --- ## Phenotypes for(i in 1:n) { y[i] ~ dnorm(mu[i], tau2) mu[i] } ## Location prior a ~ dnorm(0, 1.0E-6) b ~ dnorm(0, 1.0E-6) ## Variance prior lsigma ~ dunif(-10, 10) ## sigma between 4.539993e-05 and 2.202647e+04 tau2 sigma } modelD { ## --- Additive + Dominance model (genotype as a factor) --- ## Phenotypes for(i in 1:n) { y[i] ~ dnorm(mu[i], tau2) mu[i] } ## Location prior g[1] ~ dnorm(0, 1.0E-6) g[2] ~ dnorm(0, 1.0E-6) g[3] ~ dnorm(0, 1.0E-6) ## Variance prior lsigma ~ dunif(-10, 10) ## sigma between 4.539993e-05 and 2.202647e+04 tau2 sigma } tmpAA tmpAD tmpDA tmpDD initsA initsD fitAA n.chains=3, n.burnin=5000, n.iter=10000, n.thin=10, bugs.seed=19791123, parameters=c("a", "b", "sigma"), bugs.directory="C:/Programs/BUGS/WinBUGS14_orig") ## mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff ## a 99.9 0.0 99.9 99.9 100 100 100.0 1 1500 ## b 1.0 0.0 0.9 1.0 1 1 1.1 1 1500 ## sigma 1.0 0.0 1.0 1.0 1 1 1.0 1 400 ## deviance 2828.3 2.4 2826.0 2826.0 2828 2829 2834.0 1 1500 ## pD = 2.9 and DIC = 2831.2 fitAD