Prof. Moysés Nascimento – 4.Article
4.Artigo: Modelagem hierárquica Bayesiana na avaliação de curvas de crescimento de suínos genotipados para o gene halotano. Ciência Rural (UFSM. Impresso), v. 44, p. 1853-1859, 2014.
1- Rotinas utilizadas no software R considerando o ajuste do modelo Logístico aos dados individuais library(R2OpenBUGS) setwd("diretório") y=read.table("endereço completo do arquivo") x=read.table("endereço completo do arquivo") Nid=nrow(y) # número de animais do grupo N=ncol(y) # número de observações por animal linemodel<-function() #especificando o modelo { for (i in 1:Nid) { for (j in 1:N) { Y[i,j]~dnorm(mu[i,j],tau) mu[i,j]<-a[i]*exp(-b[i]*exp(-c[i]*X[i,j])) # modelo logístico } a[i]~dnorm(mua,taua) #priori para ai b[i]~dnorm(mub,taub) #prior para bi k[i]~dnorm(muk,tauk) #prior para ki } mua~dnorm(mu(mua),tau(mua)) # priori para mua mub~dnorm(mu(mub),tau(mub)) # priori para mub muk~dnorm(mu(muk),tau(muk)) # priori para muk taua~dgamma(mu(taua),tau(taua)) #priori para taua sigma2a<-1/sqrt(taua) taubb~dgamma(mu(taub),tau(taub)) #priori para taub sigma2a<-1/sqrt(taub) tauak~dgamma(mu(tauk),tau(tauk)) #priori para tauk sigma2a<-1/sqrt(tauk) tau~dgamma(mu(tau),tau(tau)) #priori para tau sigma2a<-1/sqrt(tau) } # Preparando os dados em formato de lista linedata<-list(Y=as.matrix(y),X=as.matrix(x),Nid = número de animais do grupo, N = número de observações por animal) lineinits=function(){list(a=rep(a1,Nid),b=rep(b1,Nid),k=rep(k1,Nid),mua= mua1, mub= mub1, muk= muk1, taua= taua1, taub= taub1,tauk= tauk1,tau= tau1)} parameters=c("a","b","c","mua","mub","muk","sigma2a","sigma2b","sigma2c","sigma2") lineout<-bugs(linedata,lineinits,parameters,linemodel, n.chains=1, n.iter=n1, n.burnin=n2, n.thin=n3, DIC=TRUE) lineout$summary #parameter estimates 2-Rotinas utilizadas no software R considerando o ajuste do modelo Logístico aos dados individuais de peso-idade dos suínos via metodologia Frequentista. #lendo dados longitudinais (peso x idade) setwd("diretório") dados=read.table("endereço completo do arquivo ",h=T) #ajuste individual das curvas de crescimento step1=by(dados, dados$id, function(x) nls(y ~ a/(1+b*exp(-k*x)), start = list(a=a1,b=b1, k=k1), data=x)) sapply(step1, coef) #arquivando as estimativas dos parâmetros a=as.matrix(sapply(step1, coef)[1,]) b=as.matrix(sapply(step1, coef)[2,]) k=as.matrix(sapply(step1, coef)[3,]) genotipos=read.table"endereço completo do arquivo ",h=T) #lendo o arquivo de genótipos M=as.matrix(genotipos[,-(1)]) #definindo matriz de marcadores trat=paste(M[,1],M[,2],sep="") table(trat) head(a) head(b) head(k) a_anova=lm(a~factor(trat)) anova(a_anova) b_anova=lm(b~factor(trat)) anova(b_anova) k_anova=lm(k~factor(trat)) anova(k_anova) |