Não foi possível enviar o arquivo. Será algum problema com as permissões?
Diferenças

Diferenças

Aqui você vê as diferenças entre duas revisões dessa página.

Link para esta página de comparações

Ambos lados da revisão anterior Revisão anterior
Próxima revisão
Revisão anterior
pessoais:eder [2011/09/25 20:44]
eder [section 2]
pessoais:eder [2018/07/04 19:10] (atual)
eder [section 5]
Linha 14: Linha 14:
   * {{:​pessoais:​inlarrblup.r|GWS}} Seleção Genômica Ampla Via ML REML INLA   * {{:​pessoais:​inlarrblup.r|GWS}} Seleção Genômica Ampla Via ML REML INLA
   * {{:​pessoais:​reml_inla.r|Script}} Modelo seleção Genótipo ambiente via REML ML INLA   * {{:​pessoais:​reml_inla.r|Script}} Modelo seleção Genótipo ambiente via REML ML INLA
 +  * {{:​pessoais:​linearregression.rnw|Script}} Regressão Linear - inferência via Mínimos quadrados, ML, REML, Gibbs, Metropolis, INLA, dclone ... (Em construção)
 +  * {{:​pessoais:​rjmcmc.r|RJMCMC}} Reversible Jump MCMC Regressão Linear ​
 +===== Artigos de Interesse ===== 
 +  * {{http://​www.sciencedirect.com/​science/​article/​pii/​S0022030278835905|Simulation of Examine Distributions of Estimators of Variances and Ratios of Variances}}
 +  * {{http://​www.jstor.org/​stable/​3001853|Estimation of Variance and Covariance Components}}
 +  * {{http://​www.sciencedirect.com/​science/​article/​pii/​S0022030291786013|C. R. Henderson: Contributions to Predicting Genetic Merit}}
 +  * {{http://​www.sciencedirect.com/​science/​article/​pii/​S002203027584776X|Rapid Method for Computing the Inverse of a Relationship Matrix}}
 +  * {{http://​www.jstor.org/​stable/​2527669?​seq=18|The Estimation of Environmental and Genetic Trends from Records Subject to Culling}}
 +  * {{http://​download.journals.elsevierhealth.com/​pdfs/​journals/​0022-0302/​PIIS002203027584776X.pdf|Rapid Method for Computing the Inverse of a Relationship Matrix}}
 +  * {{http://​www.jstor.org/​pss/​2529339|A simple method for computing the inverse of a numerator relationship matrix used in prediction of breeding values}}
 +  * {{http://​www.jstor.org/​stable/​2529430?&​Search=yes&​searchText=%22C.+R.+Henderson%22&​list=hide&​searchUri=%2Faction%2FdoBasicSearch%3FQuery%3Dau%253A%2522C.%2BR.%2BHenderson%2522%26wc%3Don&​prevSearch=&​item=3&​ttl=373&​returnArticleService=showFullText|Best Linear Unbiased Estimation and Prediction under a Selection Model}}
 +  * {{http://​www.jstor.org/​stable/​2530609?&​Search=yes&​searchText=%22C.+R.+Henderson%22&​list=hide&​searchUri=%2Faction%2FdoBasicSearch%3FQuery%3Dau%253A%2522C.%2BR.%2BHenderson%2522%26wc%3Don&​prevSearch=&​item=5&​ttl=373&​returnArticleService=showFullText|Variance-Covariance Matrix of Estimators of Variances in Unweighted Means ANOVA}}
 +  * {{http://​www.jstor.org/​stable/​3001853?&​Search=yes&​searchText=%22C.+R.+Henderson%22&​list=hide&​searchUri=%2Faction%2FdoBasicSearch%3FQuery%3Dau%253A%2522C.%2BR.%2BHenderson%2522%26wc%3Don&​prevSearch=&​item=2&​ttl=373&​returnArticleService=showFullText| Estimation of Variance and Covariance Components}}
 +
 +
 ===== Disciplinas 2011/1 =====  ===== Disciplinas 2011/1 ===== 
   * [[http://​www.leg.ufpr.br/​doku.php/​disciplinas:​ce210-2010-02|CE-210:​ Inferência estatística II]]   * [[http://​www.leg.ufpr.br/​doku.php/​disciplinas:​ce210-2010-02|CE-210:​ Inferência estatística II]]
Linha 23: Linha 38:
    * [[http://​www.leg.ufpr.br/​doku.php/​pessoais:​eder:​planejamentofito|Planejamento de experimento PG Produção Vegetal UFPR]]    * [[http://​www.leg.ufpr.br/​doku.php/​pessoais:​eder:​planejamentofito|Planejamento de experimento PG Produção Vegetal UFPR]]
    * [[http://​www.leg.ufpr.br/​doku.php/​pessoais:​eder:​exptempo| Análise de Experimentos de longa duração]] II Reunião Paranaense Ciência do Solo    * [[http://​www.leg.ufpr.br/​doku.php/​pessoais:​eder:​exptempo| Análise de Experimentos de longa duração]] II Reunião Paranaense Ciência do Solo
-===== Códigos ===== +   * [[http://​www.leg.ufpr.br/​doku.php/​pessoais:​eder:​runicentro| Curso de Sofware R Analise Experimentos - UNICENTRO]] 
 +===== Códigos ​(Em construção) ​===== 
 <code R> <code R>
-###-----------------------------------------------------------------###​ +##​------------------------------------------------------------------###​
-### Agulha de buffon +
-buffon <- function(n,​l=1,​a=1){ +
-  if(a<​l){cat('​Erro:​ a < l, deve ser a > l\n'​)} +
-  if(a>​=l){ +
-  theta <- runif(n,​0,​pi) +
-  dist <- runif(n,​0,​a/​2) +
-  inter <- sum(dist <= l/​2*sin(theta)) +
-  phi_est <- round((n/​inter)*(2*l/​a),​12) +
-  cat('​Número Simulação',​n,'​phi_estimado',​phi_est,'​Erro',​round(pi-phi_est,​12),'​\n'​) +
-  return(c(n,​phi_est)) +
-}} +
- +
-n <- seq(10000,​1000000,​by=20000) +
-res <- matrix(NA,​ncol=2,​nrow=length(n)) +
-con <- 1 +
-for (i in n){ +
-  res[con,] <- buffon(i) +
-  con <- con+1 +
-+
- +
-plot(res,​type='​l',​ylab=expression(pi),​xlab='​Simulações'​) +
-abline(h=pi,​col='​red'​) +
-###​-----------------------------------------------------------------###​ +
-### MOnte carlo +
-## Calcula a área via simulação de monte carlo +
-## args: r= raio, s vetor com numero de simulação,​ plotS plotar a simulação +
-MCcirculo<​-function(r,​s,​plotS=TRUE){ +
-ns<​-area<​-s +
-r<-r +
-con <- 1 +
-for (j in ns) { +
-#pontos aleatorios +
- x<​-runif(j,​ min=-r, max=r) +
- y<​-runif(j,​ min=-r, max=r) +
- ponto<​-cbind(x,​y) +
-  cont <- sum(apply(ponto,​1,​function(x){sqrt(sum(x^2))})<​r) +
-#plotando Simulação +
-  if(plotS==TRUE){ +
- plot(x,​y,​col="​red",​type="​p",​asp=1,​lwd=1,​xlim=c(-r,​r),​ylim=c(-r,​r),​ main="​Simulação Monte Carlo",​sub=j) +
- ang <- seq(0, 2*pi, length = 100) +
- xx <- r * cos(ang);yy <- r * sin(ang) +
- polygon(xx,​ yy,border = "dark blue",​lwd=2) +
-  }   +
-#Calculo de Area +
- area[con]<​-(cont/​j)*(r^2)*4 +
-  cat(paste(round(area[con],​6),​j,'​\n'​)) +
-  con <- con+1 +
-+
- plot(ns,​area,​main="​Simulação Monte Carlo",​xlab='​Número da amostra',​ylab='​Area'​) +
-  abline(h=pi*r^2,​col='​red',​lwd=2) +
-   +
-+
-MCcirculo(1,​seq(5,​5000,​by=1000),​plotS=FALSE) +
-###​-----------------------------------------------------------------###​ +
-### Inversão de Probabilidade +
-### OBJ: gerar x~exp transformando de uma uniforme +
-NS <- 10000 +
-lam <- 0.5 +
-#​f(x)=exp(lam) F(x)=1-exp(-lam*x),​ logo: F^-1(x)= -lam^-1*log(1-x) +
-Gexp <- function(x,​lam){-(log(1-U))/​lam} +
- +
-U <- runif(NS) +
-X <- Gexp(U,​lam) +
-Y <- rexp(NS,​lam) +
- +
-par(mfrow=c(1,​3)) +
-hist(U,​freq=FALSE,​main='​Uniforme',​col='​lightblue'​) +
-lines(density(U),​col='​red',​lwd=2) +
- +
-hist(X,​freq=FALSE,​main='​Expoencial via uniforme',​col='​lightblue'​) +
-lines(density(X),​col='​red',​lwd=2) +
-lines(curve(dexp(x,​lam),​min(X),​max(X),​add=TRUE),​col='​blue',​lwd=2) +
- +
-hist(Y,​freq=FALSE,​main='​Expoencial do R',​col='​lightblue'​) +
-lines(density(Y),​col='​red',​lwd=2) +
-lines(curve(dexp(x,​lam),​min(Y),​max(Y),​add=TRUE),​col='​blue',​lwd=2) +
-###​-----------------------------------------------------------------###​ +
-### Metodos de integração numerica +
-#Função +
-f <- function(x){exp(-x^2)} +
-a <- -3 +
-b <- 3 +
-# integrar de -3,3 +
-x <- seq(a,​b,​l=100) +
-plot(x,​f(x),​type='​l',​ylim=c(0,​1)) +
-# Integração nativa do R - Gauss–Kronrod quadrature +
-integrate(f,​a,​b) +
-###Simpson 1/3 - INtervalos par, igualmente espaçados +
-n <- 1200 +
-xi <- seq(a,​b,​l=n+1) +
-i <- seq(2,​n,​by=2) +
-j <- seq(3,​n-1,​by=2) +
-((b-a)/​n/​3)*(f(a)+4*sum(f(xi[i]))+2*sum(f(xi[j]))+f(b)) +
-###Simpson 3/8 - Intervalos divisiveis por 3 +
-n <- 1200 +
-xi <- seq(a,​b,​l=n+1) +
-i <- seq(2,​n,​by=3) +
-j <- seq(4,​n-2,​by=3) +
-((3*(b-a)/​n)/​8)*(f(a)+3*sum(f(xi[i])+f(xi[i+1]))+2*sum(f(xi[j]))+f(b)) +
-### Quadratura gausiana 3º Ordem +
-w <- c(0.555555,​0.888888,​0.555555) +
-xi <- c(-0.77459667,​0,​0.77459667) +
-(b-a)/​2*sum(f((b-a)/​2*xi+(a+b)/​2)*w) +
-### Quadratura gausiana 4º Ordem +
-w <- c(0.3478548,​0.6521452,​0.6521452,​0.3478548) +
-xi <- c(-0.86113631,​-0.33998104,​0.33998104,​0.86113631) +
-(b-a)/​2*sum(f((b-a)/​2*xi+(a+b)/​2)*w) +
-### Quadratura gausiana 6º Ordem +
-w <- c(0.1713245,​0.3607616,​0.4679139,​0.4679139,​0.3607616,​0.1713245) +
-xi <- c(-0.933246951,​-0.66120938,​-0.23861919,​0.23861919,​0.66120938,​0.933246951) +
-(b-a)/​2*sum(f((b-a)/​2*xi+(a+b)/​2)*w) +
-###Monte Carlo +
-n <- 10000 +
-xi <- runif(n,​a,​b) +
-Ls <- max(f(seq(a,​b,​l=100))) +
-Li <- 0 +
-yi <- runif(n,​Li,​Ls) +
-sum(f(xi)>​=yi)/​n*((b-a)*(Ls-Li)) +
-points(xi,​yi) +
-###​Laplace +
-#f' <- -2*x*exp(-x^2) +
-D2f  <- function(x){(4*x^2-2)*exp(-x^2)} +
-D2f(0) +
-((2*pi)/​((-D2f(0))))^0.5*f(0) +
-##​Avaliando +
-x <- seq(a,​b,​l=100) +
-plot(x,​f(x),​type='​l',​ylim=c(0,​2)) +
-lines(x,​((2*pi)/​((-D2f(0))))^0.5*f(x),​col="​red"​) +
-###​------------------------------------------------------------###​ +
-###​------------------------------------------------------------###​ +
-### Solução analitica, númerica e por simulação do modelo +
-# X ~ B(n,p) +
-# p ~ Beta(alfa,​beta) +
-###​------------------------------------------------------------###​ +
-###​------------------------------------------------------------###​ +
-require(sfsmisc) +
-require(latticeExtra) +
-require(MASS) +
-#​browseURL('​http://​cs.illinois.edu/​class/​sp10/​cs598jhm/​Slides/​Lecture02HO.pdf'​) +
- +
-###​------------------------------------------------------------###​ +
-###​------------------------------------------------------------###​ +
-### grid de p +
-p <- seq(0,​0.99999,​by=0.001) +
-### Priori +
-alfa <- 1 +
-beta <- 1 +
-p.priori <- dbeta(p,​alfa,​beta) +
-### Verossimilhança +
-n <- 1000 +
-x <- rbinom(1,​n,​0.3) +
-vero <- function(p,​n,​x){exp(sum(dbinom(x,​n,​p,​log=TRUE)))} +
-p.vero <- apply(matrix(p),​1,​vero,​n=n,​x=x) +
-###​------------------------------------------------------------###​ +
-###​------------------------------------------------------------###​ +
-### Solução analitica +
-### Posteriori +
-p.posteA <- dbeta(p,​alfa+sum(x),​beta+sum(n-x)) +
-### Plotando +
-doubleYScale(xyplot(p.priori + p.posteA ~ p, foo, type = "​l",​lwd=3),​  +
-             ​xyplot(p.vero ~ p, foo, type = "​l",​lwd=2,​lty=2),​ +
-             ​style1 = 0, style2 = 3, add.ylab2 = TRUE, +
-             text = c("​Priori",​ "​Posteriori",​ "​Verossimilhança"​),​ columns = 3)  +
-### confirmando se a posteriori é uma fdp              +
-integrate.xy(p,​p.posteA) ​             +
-###​------------------------------------------------------------###​ +
-###​------------------------------------------------------------###​ +
-### INtegração númerica para normalização +
-### posteriori +
-p.posteN <- (p.priori*p.vero)/​(integrate.xy(p,​p.priori*p.vero)) +
-### Plotando +
-doubleYScale(xyplot(p.priori + p.posteN ~ p, foo, type = "​l",​lwd=2),​  +
-             ​xyplot(p.vero ~ p, foo, type = "​l",​lwd=2,​lty=2),​ +
-             ​style1 = 0, style2 = 3, add.ylab2 = TRUE, +
-             text = c("​Priori",​ "​Posteriori",​ "​Verossimilhança"​),​ columns = 3) +
-### confirmando se a posteriori é uma fdp              +
-integrate.xy(p,​p.posteN) +
-###​------------------------------------------------------------###​ +
-###​------------------------------------------------------------###​ +
-### Amostragem da posteriori +
-ns <- 100000 +
-theta_chapeu <- sum(x)/​(n*length(x)) +
-theta_i <- rbeta(ns,​alfa,​beta) +
-    u_i <- runif(ns,​0,​1) +
-crite <- u_i <= ((dbeta(theta_i,​alfa,​beta)*apply(matrix(theta_i),​1,​vero,​n=n,​x=x))/​ +
-                 ​(dbeta(theta_chapeu,​alfa,​beta)*vero(theta_chapeu,​n=n,​x=x))) +
-a.posteriori <- theta_i[crite] ​     +
-mean(a.posteriori,​na.rm=TRUE) +
-### Taxa Aceitação +
-sum(crite)/​ns +
-###​------------------------------------------------------------###​ +
-###​------------------------------------------------------------###​ +
-### Comparando os resultados +
-hist(a.posteriori,​prob=TRUE) +
-rug(a.posteriori) +
-lines(density(a.posteriori)) +
-lines(p,​p.posteA,​col='​red',​lwd=3) +
-lines(p,​p.posteN,​col='​blue',​lty=2) +
-legend('​topleft',​c('​Amostragem','​Analitico','​Númerica'​),​lty=c(1,​1,​2),​col=c('​black','​red','​blue'​)) +
- +
-### Intervalos via verosimilhança aproximado +
-theta_chapeu+c(-1,​1)*1.96*sqrt((theta_chapeu*(1-theta_chapeu))/​n) +
-### IC amostragem +
-quantile(a.posteriori,​c(0.025,​0.975)) +
-### Analitico da conjugada +
-qbeta(c(0.025,​0.975),​alfa+sum(x),​beta+sum(n-x)) +
-###​------------------------------------------------------------###​ +
-##------------------------------------------------------------###​+
 ###​-----------------------------------------------------------------###​ ###​-----------------------------------------------------------------###​
 ### Regressão Beta ### Regressão Beta
Linha 277: Linha 85:
 plot(grid.phi,​2*(max(opt)-opt),​type='​l'​) plot(grid.phi,​2*(max(opt)-opt),​type='​l'​)
 abline(h=3.84) abline(h=3.84)
-##############################################################################################​ 
-###​--------------------------------------------------###​ 
-###​--------------------------------------------------###​ 
-### Algumas funções para analise utilizando MCMC 
-### topicos: 
-### 1) Modelos com JAGS via R (Regressão linear, Logistica) 
-### 2) dclone juntamente com JAGS 
-### 3) MCMCsamp em modelos mistos 
-### 4) MCMCglmm 
-### 5) MCMCpack 
-###​--------------------------------------------------###​ 
-###​--------------------------------------------------###​ 
-rm(list=ls()) 
-### JAGS 
-require(runjags) 
-###​--------------------------------------------------###​ 
-#### Exemplo JAGS (função run-jags) 
-## Exemplo de regressão linear simples 
-# Simulação de dados 
-X <- 1:50 
-Y <- rnorm(length(X),​ 2*X + 10, 6) 
-plot(X,Y) 
-###​--------------------------------------------------###​ 
-### Inferencia por minimos quadrados 
-m0 <- lm(Y~X) 
-summary(m0) 
-###​--------------------------------------------------###​ 
-### Modelo y = a*x+b 
-### Ajustando com o JAGS 
-### Priori para a ~ dnorm(0,​0.001) 
-### Priori para b ~ dnorm(0,​0.001) 
-### Priori para sigma ~ exp(1) 
-### Escrevendo o modelo 
-model <- "model {  
-for(i in 1 : N){  
-   Y[i] ~ dnorm(y.est[i],​sigma); ​ 
-   ​y.est[i] <- (a * X[i]) + b; 
- 
-a ~ dnorm(0,​0.001); ​ 
-b ~ dnorm(0,​0.001);​ 
-sigma ~ dexp(1); ​ 
-}" 
- 
-### Data e valor inicial 
-data <- dump.format(list(X=X,​ Y=Y, N=length(X))) 
-inits0 <- dump.format(list(m=1,​ c=10, sigma=1)) 
-inits1 <- dump.format(list(m=0,​ c=0, sigma=1)) 
-# Run the model  
-m1 <- run.jags(model=model,​ monitor=c("​a","​b","​sigma"​), ​ 
-data=data, inits=c(inits0,​inits1),​n.chains=3,​ plots = TRUE) 
-### informações do objeto 
-names(m1) 
-### Verificando a cadeia 
-plot(m1$mcmc[[1]]) 
-### Intervalos de credibilidade 
-m1$HPD 
-### Sumario 
-m1$summary 
-###​--------------------------------------------------###​ 
-###​--------------------------------------------------###​ 
-### Regressão logistica (Livro Introdução Analise Bayesiana) 
-## pi (n/y) é a proporção de embriões com asas  
-## ti é a variável “tempo desde a deposição dos ovos” 
-## função de ligação logit 
- 
-t <- c(5,​6,​8,​8,​10,​11,​16,​18) 
-n <- c(34,​33,​33,​35,​30,​27,​33,​39) 
-y <- c(6,​4,​23,​18,​28,​27,​33,​39) 
- 
-dados <- cbind(y,​n-y) 
-m0 <- glm(dados~t,​family=binomial(link=logit)) 
-summary(m0) 
- 
-datalist <- dump.format(list(metamorfose=y,​total=n,​tempo=t)) 
-params <- c("​beta0","​beta1"​) 
-inicial <- dump.format(list(beta0=coef(m0)[1],​ 
-                            beta1=coef(m0)[2])) 
-###​--------------------------------------------------###​ 
-modmcmc <- "​model{ 
-for(i in 1:​length(metamorfose)){ 
-    metamorfose[i] ~ dbin(p[i],​total[i]) 
-    logit(p[i]) <- beta0+beta1*tempo[i] 
-} 
-beta0 ~ dnorm(0,​0.001) 
-beta1 ~ dnorm(0,​0.001) 
-}" 
- 
-### rodando o modelo 
-modfit <- run.jags(model=modmcmc,​monitor=params,​data=datalist,​ 
-inits=inicial,​n.chains=1,​burnin=10000,​thin=3,​sample=9000,​check.conv=TRUE) 
-### Verificando a cadeia 
-plot(modfit$mcmc[[1]]) 
-### Intervalos de credibilidade 
-modfit$HPD 
-### Sumario 
-modfit$summary 
-###​--------------------------------------------------###​ 
-###​--------------------------------------------------###​ 
-### Modelos com repetição 
-x <- gl(5,4) 
-y <- rnorm(20,​10,​5) 
-X <- model.matrix(~x) 
-datalist <- dump.format(list(y=y,​X=X)) 
-params <- c("​beta1","​tau"​) 
-inicial <- dump.format(list(beta1=rep(0,​nlevels(x)),​tau=1)) 
-###​--------------------------------------------------###​ 
-modmcmc <- "​model{ 
-for(i in 1:​length(y)){ 
-    y[i] ~ dnorm(yb[i],​tau) 
-    yb[i] <- inprod(X[i,​],​ beta1) 
-} 
-tau ~ dexp(1); 
-for(j in 1:5){ 
-    beta1[j] ~ dnorm(0,​0.001) 
-} 
-}" 
-### rodando o modelo 
-modfit <- run.jags(model=modmcmc,​monitor=params,​data=datalist,​ 
-inits=inicial,​check.conv=TRUE) 
-### Verificando a cadeia 
-plot(modfit$mcmc[[1]]) 
-### Intervalos de credibilidade 
-modfit$HPD 
-### Sumario 
-modfit$summary 
-#### Minimos quadrados 
-summary(lm(y~x)) 
-###​--------------------------------------------------###​ 
-###​--------------------------------------------------###​ 
-### JAGS com Dclone 
-require(dclone) 
-## simple regression example from the JAGS manual 
-jfun <- function() { 
-for (i in 1:N) { 
-  Y[i] ~ dnorm(mu[i],​ tau) 
-  mu[i] <- alpha + beta * (x[i] - x.bar) 
-} 
-x.bar <- mean(x[]) 
-alpha ~ dnorm(0.0, 1.0E-4) 
-beta ~ dnorm(0.0, 1.0E-4) 
-sigma <- 1.0/​sqrt(tau) 
-tau ~ dgamma(1.0E-3,​ 1.0E-3) 
-} 
-## data generation 
-set.seed(1234) 
-N <- 100 
-alpha <- 1 
-beta <- -1 
-sigma <- 0.5 
-x <- runif(N) 
-linpred <- model.matrix(~x) %*% c(alpha, beta) 
-Y <- rnorm(N, mean = linpred, sd = sigma) 
-## list of data for the model 
-jdata <- list(N = N, Y = Y, x = x) 
-## what to monitor 
-jpara <- c("​alpha",​ "​beta",​ "​sigma"​) 
-## fit the model with JAGS 
-regmod <- jags.fit(jdata,​ jpara, jfun, n.chains = 3) 
-## model summary 
-summary(regmod) 
-## data cloning 
-dcdata <- dclone(jdata,​ 5, multiply = "​N"​) 
-dcmod <- jags.fit(dcdata,​ jpara, jfun, n.chains = 3) 
-summary(dcmod) 
-#### opções de computação paralela 
-?​jags.parfit 
  
-###​--------------------------------------------------###​ 
-###​--------------------------------------------------###​ 
-### Modelo misto - Amostrando a posteriorir via MCMC 
-require(lme4) 
-### Resposta normal 
-m0 <- lmer(Reaction ~ Days + (1|Subject) + (0+Days|Subject),​ sleepstudy) 
-summary(m0) 
-sampm0 <- mcmcsamp(m0,​ n = 1000) 
-HPDinterval(sampm0) 
-xyplot(sampm0) 
-qqmath(sampm0) 
-densityplot(sampm0) 
-###​--------------------------------------------------###​ 
-###​--------------------------------------------------###​ 
-browseURL("​http://​cran-r.c3sl.ufpr.br/​web/​packages/​MCMCglmm/​vignettes/​CourseNotes.pdf"​) 
-browseURL("​http://​cran-r.c3sl.ufpr.br/​web/​packages/​MCMCglmm/​vignettes/​Overview.pdf"​) 
-browseURL("​http://​citeseerx.ist.psu.edu/​viewdoc/​download?​doi=10.1.1.160.5098&​rep=rep1&​type=pdf"​) 
-require(MCMCglmm) 
-?MCMCglmm 
-data(PlodiaPO) 
-### Modelo normal 
-model1<​-MCMCglmm(PO~1,​ random=~FSfamily,​ data=PlodiaPO,​ verbose=FALSE) 
-summary(model1) 
-plot.MCMCglmm(model1) 
-### MOdelo binomial 
-model2 <- MCMCglmm(cbind(Pupated,​ Infected) ~ 1,​random=~FSfamily,​ family = "​multinomial2",​ 
-data = PlodiaR, verbose = FALSE) 
-plot.MCMCglmm(model2) 
-###​--------------------------------------------------###​ 
-###​--------------------------------------------------###​ 
-require(MCMCpack) 
-### Modelo Poisson 
-counts <- c(18,​17,​15,​20,​10,​20,​25,​13,​12) 
-outcome <- gl(3,1,9) 
-treatment <- gl(3,3) 
-posterior <- MCMCpoisson(counts ~ outcome + treatment) 
-plot(posterior) 
-summary(posterior) 
-###​--------------------------------------------------###​ 
-###​--------------------------------------------------###​ 
-##############################################################################################​ 
 </​code>​ </​code>​
  

QR Code
QR Code pessoais:eder (generated for current page)