### Seleção genomica ampla - testes #browseURL('http://www.infoteca.cnptia.embrapa.br/bitstream/doc/883425/1/Doc210.pdf') #pg 54 ###-----------------------------------------------------### rm(list=ls()) require(INLA) require(rrBLUP) require(plyr) require(lattice) ###-----------------------------------------------------### dados <- data.frame(ind=c(1:5), diametro=c(9.87,14.48,8.91,14.64,9.55), M1=c(2,1,0,1,1), M2=c(0,1,2,0,0), M3=c(0,0,0,1,0), M4=c(0,0,0,0,1), M5=c(2,1,0,1,1), M6=c(0,1,0,0,1), M7=c(0,0,2,0,0)) dados ###-----------------------------------------------------### ### Criando a MAtrix Z de efeitos aleatorios Z <- as.matrix(dados[,3:ncol(dados)]) ### Repadronizando de 0 1 2 para -1 0 1 Z <- apply(Z,2,function(x) ifelse(x==1,-1, ifelse(x!=0,1,0))) levelplot(Z) ###-----------------------------------------------------### ### Abordagem REML result.reml <- mixed.solve(y=dados$diametro, Z=Z, K=diag(ncol(Z)),method = "REML") result.reml ###-----------------------------------------------------### ### Abordagem ML REML result.ml <- mixed.solve(y=dados$diametro, Z=Z, K=diag(ncol(Z)),method = "ML") result.ml ###-----------------------------------------------------### ### Abordagem Bayesiana via INLA rr.inla.fit = inla(diametro ~ 1 + f(ind,model="z",Z=Z),data=dados, family="gaussian") summary(rr.inla.fit) 1/rr.inla.fit$summary.hyperpar[,1] rr.inla.fit$'summary.random' ###-----------------------------------------------------### ### Componentes de variancia compVAR <- c(result.reml$Vu,result.ml$Vu,(1/rr.inla.fit$summary.hyperpar[,1])[2]) names(compVAR) <- c('REML','ML','INLA') compVAR compVARe <- round(c(result.reml$Ve,result.ml$Ve,(1/rr.inla.fit$summary.hyperpar[,1])[1]),7) names(compVARe) <- c('REML','ML','INLA') compVARe ### Razao compVAR/(compVAR+compVARe) ###-----------------------------------------------------### ### Plotando os efeitos e.ml <- result.ml$u e.reml <- result.reml$u e.inla <- ldply(rr.inla.fit$'summary.random')[,3] xyplot(e.reml + e.ml + e.inla ~ 1:ncol(Z),type=c('+','o'), ylab='Efeitos',xlab='Marcadores',jitter.y=TRUE,auto.key=TRUE) ### IC INLA e.inlaIC0.025 <- ldply(rr.inla.fit$'summary.random')[,5] e.inlaIC0.975 <- ldply(rr.inla.fit$'summary.random')[,7] xyplot(e.inla + e.inlaIC0.025+e.inlaIC0.975 ~ 1:ncol(Z),type=c('+','o'), ylab='Efeitos',xlab='Marcadores',jitter.y=TRUE,auto.key=TRUE) ###-----------------------------------------------------### ###-----------------------------------------------------### ### Simulando dados Nind <- 100 Nmarc <- 500 set.seed(1) dados <- data.frame(ind=c(1:5),diametro=rnorm(Nind,50,25),matrix(sample(-1:1,Nind*Nmarc,replace =TRUE),nrow=Nind,ncol=Nmarc)) dados ###-----------------------------------------------------### ### Criando a MAtrix Z de efeitos aleatorios Z <- as.matrix(dados[,3:ncol(dados)]) levelplot(Z,xlab='',ylab='') ###-----------------------------------------------------### ### Abordagem REML result.reml <- mixed.solve(y=dados$diametro, Z=Z, K=diag(ncol(Z)),method = "REML") result.reml ###-----------------------------------------------------### ### Abordagem ML REML result.ml <- mixed.solve(y=dados$diametro, Z=Z, K=diag(ncol(Z)),method = "ML",SE=TRUE) result.ml ###-----------------------------------------------------### ### Abordagem Bayesiana via INLA rr.inla.fit = inla(diametro ~ 1 + f(ind,model="z",Z=Z),data=dados, family="gaussian") summary(rr.inla.fit) 1/rr.inla.fit$summary.hyperpar[,1] rr.inla.fit$'summary.random' ###-----------------------------------------------------### ### Componentes de variancia compVAR <- c(result.reml$Vu,result.ml$Vu,(1/rr.inla.fit$summary.hyperpar[,1])[2]) names(compVAR) <- c('REML','ML','INLA') compVAR compVARe <- round(c(result.reml$Ve,result.ml$Ve,(1/rr.inla.fit$summary.hyperpar[,1])[1]),7) names(compVARe) <- c('REML','ML','INLA') compVARe ### Razao round(compVAR/(compVAR+compVARe),6) ###-----------------------------------------------------### ### Plotando os efeitos e.ml <- result.ml$u e.reml <- result.reml$u e.inla <- ldply(rr.inla.fit$'summary.random')[,3] xyplot(e.reml + e.ml + e.inla ~ 1:ncol(Z),type='l', ylab='Efeitos',xlab='Marcadores',jitter.y=TRUE,auto.key=TRUE) ### IC INLA e.inlaIC0.025 <- ldply(rr.inla.fit$'summary.random')[,5] e.inlaIC0.975 <- ldply(rr.inla.fit$'summary.random')[,7] xyplot(e.inla+e.inlaIC0.025+e.inlaIC0.975 ~ 1:ncol(Z),type='l', ylab='Efeitos',xlab='Marcadores',jitter.y=TRUE,auto.key=TRUE) ###-----------------------------------------------------###