• Graduado em Engenharia Agronômica - Eng Agronômica - UFPR

• e-mail: ederdbs@gmail.com / eder@leg.ufpr.br

• Currículo LATTES
• Arquivos

• Estatística Experimental
• Estatística Espacial
• GEM² Grupo de estudos em modelos mistos
• GWS Seleção Genômica Ampla Via ML REML INLA
• Script Modelo seleção Genótipo ambiente via REML ML INLA
• Script Regressão Linear - inferência via Mínimos quadrados, ML, REML, Gibbs, Metropolis, INLA, dclone … (Em construção)

• CE-210: Inferência estatística II
• CE-718: Métodos Computacionalmente Intensivos

• Estatística Experimental com Software R
• Planejamento de experimento PG Produção Vegetal UFPR
• Análise de Experimentos de longa duração II Reunião Paranaense Ciência do Solo

###-----------------------------------------------------------------###
### Reversible jump MCMC
### Modelo 1 y ~ N(b0+b1*x,sigma)
### Modelo 1 y ~ N(b0+b1*x+b2*x²,sigma)
#browseURL('http://www.icmc.usp.br/~ehlers/bayes/cap4.pdf')
# pg 76
require(MASS)#mvnorm()
require(MCMCpack)#rinvgamma() dinvgamma()
require(coda)#as.mcmc
rm(list=ls())
### Ajustar Prioris
### conferir jacobiano
 
rj.modelo <- function(y,x,b0,b1,sigma,b01,b11,b21,sigma1,model,mu=mu,sd=sd,mu0=mu0,
                      mu02=mu02,V0=V0,V02=V02,v0=v0,tau0=tau0,v02=v02,tau02=tau02){          
    if (model == 1){
      u <- rnorm(1, mu,sd)
      b0_n <- b0 
      b1_n <- b1
      sigma_n <- sigma
      b01_n <- b0 * u  
      b11_n <-  b1 * u 
      b21_n <- u
      sigma1_n <- sigma * (u^2)
      }
    if (model == 2){
      u <- b21
      b0_n <- b01 / u
      b1_n <- b11 / u
      sigma_n <-  sigma1 / (u^2)
      b01_n <- b01
      b11_n <-  b11
      b21_n <- b21
      sigma1_n <- sigma1        
      }
    num <-  (sum(dnorm(y,b0_n+b1_n*x,sigma_n,log=TRUE))#+
             #sum(dnorm(b0_n,mu0[1],V0[1,1],log=TRUE))+
             #sum(dnorm(b1_n,mu0[1],V0[2,2],log=TRUE))+
             #sum(log(dinvgamma(sigma_n,v0,tau0)))
             ) * u^4
    den <-  (sum(dnorm(y,b01_n+b11_n*x+b21_n*x^2,sigma1_n,log=TRUE))#+
            # sum(dnorm(b01_n,mu02[1],V02[1,1],log=TRUE))+
            # sum(dnorm(b11_n,mu02[2],V02[2,2],log=TRUE))+
            # sum(dnorm(b21_n,mu02[3],V02[3,3],log=TRUE))+
            # sum(log(dinvgamma(sigma1_n,v02,tau02)))
             ) * dnorm(u,0,2)
      u = runif(1, 0, 1)
      if (model == 1) {
        aceita = min(1, num/den)
          if (u < aceita) {
          model = 2
          b0 <- b0_n
          b1 <- b1_n
          sigma <- sigma_n
          }
      }
      if (model == 2){
          aceita = min(1, den/num)
            if (u < aceita) {
              model = 1
              b01 <- b01_n
              b11 <- b11_n
              b21 <- b21_n
              sigma1 <- sigma1_n            
              }
          }
      if (model == 1){return(list(model = model,b0=b0,b1=b1,sigma=sigma))}
      if (model == 2){return(list(model = model,b01=b01,b11=b11,b21=b21,sigma1=sigma1))}
  }
 
 
 
rjmcmc <- function(nI, x,y,burnIN,mu=mu,sd=sd) {
  chain = matrix(NA, nrow = nI, ncol = 8)
  nv <- c(0,0)
  chain[1,1:8] = c(1)
  model = 1
  n <- length(y)
  ###----------------------------------------------------------###
  ###MOdel 1
  X <- model.matrix(~x)
  k<-ncol(X)
  #beta
  mu0<-rep(0,k)
  V0<-100*diag(k)
  #sigma2
  v0<-3
  tau0<-100
  #Valores iniciais
  chain[1,3] <-sig2draw<- 3
  invV0 <- solve(V0)  
  XtX <- crossprod(X,X)
  Xty <- crossprod(X,y)
  invV0_mu0 <- invV0 %*% mu0
  ###----------------------------------------------------------###
  # Model 2
  X2 <- cbind(1,x,x^2)
  k2<-ncol(X2)
  #beta
  mu02<-rep(0,k2)
  V02<-100*diag(k2)
  #sigma2
  v02<-3
  tau02<-100
  #Valores iniciais
  chain[1,7] <- sig2draw2<- 3
  invV02 <- solve(V02)  
  XtX2 <- crossprod(X2,X2)
  Xty2 <- crossprod(X2,y)
  invV0_mu02 <- invV02 %*% mu02  
  ###----------------------------------------------------------###
  for (i in 2:nI) {
    if (model == 1){
         # Model 1
         #beta
         invsig2draw <- 1/sig2draw
         V1<-solve(invV0+(invsig2draw) * XtX)
         mu1<-V1 %*% (invV0_mu0 + (invsig2draw)* Xty)
         chain[i,1:2]<-mvrnorm(n=1,mu1,V1)
         # sigma
         v1<-(n+2*v0)/2
         yXb <- (y-X %*% chain[i,1:2])
         tyXb <-t(yXb)
         tau1<-(0.5)*(tyXb %*% yXb+2*tau0)
         chain[i,3] <- sig2draw <- sqrt(rinvgamma(1,v1,tau1))
      }
    if (model == 2){
         # Model 2
         #beta
         invsig2draw2 <- 1/sig2draw2
         V12<-solve(invV02+(invsig2draw2) * XtX2)
         mu12<-V12 %*% (invV0_mu02 + (invsig2draw2)* Xty2)
         chain[i,4:6]<-mvrnorm(n=1,mu12,V12)
         # sigma
         v12<-(n+2*v02)/2
         yXb2 <- (y-X2 %*% chain[i,4:6])
         tyXb2 <-t(yXb2)
         tau12<-(0.5)*(tyXb2 %*% yXb2+2*tau02)
         chain[i,7] <- sig2draw2 <-  sqrt(rinvgamma(1,v12,tau12))
        }
    new <- rj.modelo(y,x,chain[i,1],chain[i,2],chain[i,3],chain[i,4],chain[i,5],chain[i,6],chain[i,7],model,mu=mu,sd=sd)
    model  <-  new$model
    if (model == 1) {
            chain[i, 1] = new$b0
            chain[i, 2] = new$b1
            chain[i, 3] = new$sigma
            nv[1] = nv[1] + 1
            }
    if (model == 2) {
            chain[i, 4] = new$b01
            chain[i, 5] = new$b11
            chain[i, 6] = new$b21
            chain[i, 7] = new$sigma1
            nv[2] = nv[2] + 1
    }
}
chain[,8] <- 1
chain[is.na(chain[,1]),8] <- 2
chain <- chain[- c(1:burnIN),]
colnames(chain) <- c('b0_1','b1_1','sigma_1','b0_2','b1_2','b2_2','sigma_2','model')
return(list(as.mcmc(na.omit(chain[,1:3])),
            as.mcmc(na.omit(chain[,4:7])),
            as.mcmc(na.omit(chain[,8]))))
}
 
x <- 1:10
y <- 10+2*x^1+rnorm(x,0,5)
plot(x,y)
res <- rjmcmc(5000,x,y,1,mu=0,sd=100)
lapply(res,summary)
plot(res[[1]])
summary(lm(y~1+I(x)))          
plot(res[[2]])
summary(lm(y~1+I(x)+I(x^2)))
plot(res[[3]])
##------------------------------------------------------------------###
###-----------------------------------------------------------------###
### Regressão Beta
### pacote oficial
require(betareg)
data("FoodExpenditure", package = "betareg")
fe_beta <- betareg(I(food/income) ~ income + persons , data = FoodExpenditure)
summary(fe_beta)
###-----------------------------------------------------------------###
### log vero da regressão beta com duas covariaveis, 
log.vero <- function(par,y,x1,x2){
        mu <- exp((par[1] + par[2] * x1 + par[3] * x2))/(1+exp((par[1] + par[2] * x1 + par[3] * x2)))##logit^-1
        ll  <- sum(dbeta(y, mu* par[4], (1-mu)*par[4],log = TRUE))
        return(ll)
}
 
###-----------------------------------------------------------------###          
opt <- optim(c(B0=-0.5,B1=-0.51,B2=0.11,phi=35),log.vero,y=FoodExpenditure$food/FoodExpenditure$income,
                                                        x1=FoodExpenditure$income,
                                                        x2=FoodExpenditure$persons,
                                                        hessian = TRUE, control=(list(fnscale=-1)))
opt
opt$par
sqrt(-diag(solve(opt$hessian)))
summary(fe_beta)
###-----------------------------------------------------------------###
log.veroP <- function(par,phi,y,x1,x2){
        mu <- exp((par[1] + par[2] * x1 + par[3] * x2))/(1+exp((par[1] + par[2] * x1 + par[3] * x2)))##logit^-1
        ll  <- sum(dbeta(y, mu* phi, (1-mu)*phi,log = TRUE))
        return(ll)
}
 
opt <- grid.phi <- seq(20,60,l=150)
con <- 1
for (i in grid.phi){
  opt[con] <- optim(c(B0=-0.5,B1=-0.51,B2=0.11),log.veroP,phi=i,y=FoodExpenditure$food/FoodExpenditure$income,
                                                        x1=FoodExpenditure$income,
                                                        x2=FoodExpenditure$persons,
                                                        hessian = TRUE, control=(list(fnscale=-1)))$value
  con <- con+1
}
 
plot(grid.phi,2*(max(opt)-opt),type='l')
abline(h=3.84)

Regressão beta