R/funcoes.R

In vitorcapdeville/modHier:

#' Comparação de média de 2 grupos
#'
#' @param y1 Observações relativas ao grupo 1.
#' @param y2 Observações relativas ao grupo 2.
#' @param mu Valor inicial para mu.
#' @param del Valor inicial para delta.
#' @param mu0 Valor da média a priori de mu.
#' @param g02 Valor da variância a priori de mu.
#' @param del0 Valor da média a priori de delta.
#' @param t02 Valor da variância a priori de delta.
#' @param s20 Valor do parâmetro 1 da priori de sigma2.
#' @param nu0 valor do parâmetro 2 da priori de sigma2.
#' @param nit Número de iterações.
#' @param seed Valor da semente aleatória.
#'
#' @return Lista contendo uma amostra da posteriori para cada parâmetro.
#' @export
#'

two_groups <- function(y1,y2,mu,del,mu0,g02,del0,t02,s20,nu0,nit = 5000,seed = NULL){
  MU <- DEL <- S2 <- NULL
  set.seed(seed)
  for(s in 1:nit){
    ##update s2
    s2<-1/rgamma (1,(nu0+n1+n2)/2,
                  (nu0*s20+sum((y1-mu-del)^2)+sum((y2-mu+del)^2))/2)
    ##
    ##update mu
    var.mu<- 1/(1/g02+(n1+n2)/s2)
    mean.mu<-var.mu*(mu0/g02+sum(y1-del)/s2+sum(y2+del)/s2)
    mu<-rnorm(1,mean.mu,sqrt(var.mu))
    ##
    ##update delta
    var.del<-1/(1/t02+(n1+n2)/s2)
    mean.del<-var.del*(del0/t02+sum(y1-mu)/s2-sum(y2-mu)/s2)
    del<-rnorm(1,mean.del,sqrt(var.del))
    ##
    ##save parameter values
    MU <- c(MU,mu);DEL<-c(DEL,del);S2<-c(S2,s2)
  }
  return(list(mu = MU, delta = DEL, s2 = S2))
}

#' Modelo hierárquico Normal
#'
#' @param Y Matriz nx2, onde n é o número de observações. A primeira coluna deve ser numerada
#' de 1 a m, onde m é o número de grupos. A segunda coluna deve conter os valores observados.
#' @param mu Valor inicial para mu. Escalar.
#' @param theta Valor inicial para theta. Vetor de tamanho m.
#' @param sigma2 Valor inicial para sigma2. Escalar.
#' @param tau2 Valor inicial para tau2. Escalar.
#' @param mu0 Valor da média a priori de mu.
#' @param g02 Valor da variância a priori de mu.
#' @param del0 Valor da média a priori de delta.
#' @param t02 Valor da variância a priori de delta.
#' @param s20 Valor do parâmetro 1 da priori de sigma2.
#' @param nu0 valor do parâmetro 2 da priori de sigma2.
#' @param nit Número de iterações.
#' @param seed Valor da semente aleatória.
#'
#' @return Lista contendo uma amostra da posteriori para cada parâmetro.
#' @export
#'

hier_norm <- function(Y,mu,theta,sigma2,tau2,mu0,g02,del0,t02,s20,nu0,nit = 5000,seed = NULL){
  m<-length(unique(Y[,1]))
  n <- rep(NA,m)
  for(j in 1:m){
    n[j]<-sum(Y[,1]==j)
  }
  set.seed(seed)
  THETA<-matrix(nrow=nit,ncol=m)
  SMT<-matrix(nrow=nit,ncol=3)
  ###
  ###MCMC algorithm
  for(s in 1:nit){
    # sample new values of the thetas
    for(j in 1:m){
      vtheta<-1/(n[j]/sigma2+1/tau2)
      etheta<-vtheta*(ybar[j]*n[j]/sigma2+mu/tau2)
      theta[j]<-rnorm(1,etheta,sqrt(vtheta))
    }
    #sample new value o f sigma2
    nun<-nu0+sum(n)
    ss<-nu0*s20
    for(j in 1:m){ss<-ss+sum((Y[Y[,1]==j,2]-theta[j])^2)}
    sigma2<-1/rgamma(1,nun/2,ss/2)
    #sample a new value of mu
    vmu<- 1/(m/tau2+1/g02)
    emu<-vmu*(m*mean(theta)/tau2+mu0/g02)
    mu<-rnorm(1,emu,sqrt(vmu))
    # sample a new value of tau2
    etam<-eta0+m
    ss<-eta0*t02+sum((theta-mu)^2)
    tau2<-1/rgamma(1,etam/2,ss/2)
    #store results
    THETA[s,]<-theta
    SMT[s,]<-c(sigma2,mu,tau2)
  } ###
  return(list(theta = THETA, sigma2 = SMT[,1], mu = SMT[,2], tau2 = SMT[,3]))
}

#' Modelo efeitos mistos logístico
#'
#' @param y vetor com a quantidade de sucessos para cada indivíduo.
#' @param ni vetor com a quantidade de tentativas para cada indivíduo.
#' @param x matriz nxp com as covariáveis.
#' @param rb desvio padrão da proposta para beta.
#' @param rg vetor de desvio padrão para a prosposta para gama. Caso não seja fornecido,
#'  o desvio padrão é 1 para todos os gamas
#' @param beta.init Valores iniciais para beta. Os valores inciais para gama são fixos em 0.
#' @param nit Número de iterações.
#' @param seed Valor da semente aleatória.
#'
#' @return Lista contendo uma amostra da posteriori para cada parâmetro.
#' @export
#'

ef_mist_logit <- function(y,ni,x,rb= 0.2,rg=NULL,beta.init = NULL,nit = 5000,seed = NULL){
  n = length(y)
  set.seed(seed)
  if(is.null(beta.init)){
    beta.init = glm(cbind(y,ni-y)~x[,-1],family=binomial("logit"))$coef
  }
  if(is.null(rg)){
    rg = rep(1,n)
  }
  p = ncol(x)
  M = nit
  a2 = b2 = 0.1
  phi2 = 1;
  gam = matrix(NA,M,n); gam[1,] = 0;
  bet =  matrix(NA,M,p); bet[1,] = beta.init
  ac1 = 0;ac2 = rep(0,n);
  for (k in 2:M){

    ### Bloco phi2
    phi2[k] = rgamma(1,(n+a2)/2,(sum(gam[k-1,]^2)+b2)/2)

    ### Bloco beta
    bet.prop = bet[k-1,]+rb*rnorm(p)

    num.bet = lvero(bet.prop,gam[k-1,],x,y,ni)+sum(dnorm(bet.prop,0,10,log=T))
    den.bet = lvero(as.vector(bet[k-1,]),as.vector(gam[k-1,]),x,y,ni)+sum(dnorm(bet[k-1,],0,10,log=T))
    r.bet = num.bet - den.bet

    u = runif(1)
    if (log(u)<r.bet){bet[k,] = bet.prop; ac1 = ac1 + 1;} else{bet[k,] = bet[k-1,];}

    ### Bloco gamma
    for (i in 1:n){
      gam.prop = gam[k-1,i] + rg[i]*rnorm(1);

      num.gam = lveroi(bet[k,],gam.prop,x,y,ni,i)+dnorm(gam.prop,0,sqrt(1/phi2[k]),log=T)
      den.gam = lveroi(bet[k,],gam[k-1,i],x,y,ni,i)+dnorm(gam[k-1,i],0,sqrt(1/phi2[k]),log=T)
      r.gam = num.gam - den.gam

      u = runif(1)
      if (log(u)<r.gam){gam[k,i] = gam.prop; ac2[i] = ac2[i] + 1;} else{gam[k,i] = gam[k-1,i];}
    }
    if(k%%100 == 0){print(k)}
  }
  return(list(phi2 = phi2,beta = bet,gama=gam))
}

#' verossimilhança para a observação i

lveroi = function(bet,gami,x,y,ni,i){
  aux = exp(x[i,]%*%bet+gami)
  aux.pi = aux/(1+aux)
  result = y[i]*log(aux.pi)+(ni[i]-y[i])*log(1-aux.pi)
  #result = dbinom(y,ni,prob=aux.pi,log=T)
  result
}

#' verossimilhança para o vetor de observações

lvero = function(bet,gam,x,y,ni){
  aux = exp(x%*%bet+gam)
  aux.pi = aux/(1+aux)
  result = y*log(aux.pi)+(ni-y)*log(1-aux.pi)
  #result = dbinom(y,ni,prob=aux.pi,log=T)
  sum(result)
}

#' Modelo efeitos mistos
#'
#' \deqn{y_ij ~ N(mu_ij, 1/phi1)}
#' \deqn{mu_ij = alpha + beta*x[ij] + gamma[i]}
#' \deqn{phi1 ~ gamma(a1/2,b1/2)}
#' \deqn{alpha ~ N(0,c1)}
#' \deqn{beta ~ N(0,c2)}
#' \deqn{gamma[i] | phi2 ~ N(0,1/phi2)}
#' \deqn{phi2 ~ gamma(a2/2,b2/2)}
#'
#' @param Y matriz com os valores da variável resposta em cada momento no tempo.
#' As linhas representam os momentos no tempo e as colunas representam os indivíduos.
#' @param x matriz com a mesma estrutura que a matriz Y.
#' @param nit Número de iterações.
#' @param seed Valor da semente aleatória.
#'
#' @return Lista contendo uma amostra da posteriori para cada parâmetro.
#' @export
#'

ef_mist = function(Y, x, a1 = 0.01, b1 = 0.01, c1 = 100, c2 = 100, a2 = 0.01, b2 = 0.01, init = c(0.5,0.5,0,0),nit = 5000,seed=NULL){
  set.seed(seed)
  n = ncol(Y); J = nrow(Y);
  nobs = 0;
  for (i in 1:n){
    aux = (Y[,i]!="NA");
    nobs[i] = sum(aux,na.rm=T)
  }
  N = sum(nobs)
  xnew = NA*x
  for (i in 1:n){
    xnew[1:nobs[i],i] = x[1:nobs[i],i]
  }
  sx = sum(xnew^2,na.rm=T)
  syi = apply(Y,2,sum,na.rm=T)
  sxi = apply(xnew,2,sum,na.rm=T)

  #### inicializando
  M = nit
  gam = matrix(NA,M,n)
  gam[1,] = 0
  phi1 = init[1]
  phi2 = init[2]
  bet = init[3]
  alph = init[4]


  for (k in 2:M){
    ##### Bloco phi1
    gamk = matrix(gam[k-1,],J,n,byrow=T)
    muk = gamk+alph[k-1]+bet[k-1]*xnew
    s1 = sum((Y-muk)^2,na.rm=T)
    phi1[k] = rgamma(1,(N+a1)/2,(s1+b1)/2)
    ##### Bloco beta
    v12 = 1/(phi1[k]*sx+1/c2)
    sxy = sum(x*(Y-alph[k-1]-gamk),na.rm=T)
    mu1 = v12*phi1[k]*sxy
    bet[k] = rnorm(1,mu1,sqrt(v12))
    ##### Bloco gamma
    for (i in 1:n){
      v22 = 1/(nobs[i]*phi1[k]+phi2[k-1])
      mu2 = phi1[k]*v22*(syi[i]-nobs[i]*alph[k-1]-bet[k]*sxi[i])
      gam[k,i] = rnorm(1,mu2,sqrt(v22))
    }
    ##### Bloco alpha
    gamk = matrix(gam[k,],J,n,byrow=T)
    v32 = 1/(phi1[k]*N+1/c1)
    sxy2 = sum((Y-bet[k]*xnew-gamk),na.rm=T)
    mu3 = v32*phi1[k]*sxy2
    alph[k] = rnorm(1,mu3,sqrt(v32))
    ##### Bloco phi2
    sg = sum(gam[k,]^2)
    phi2[k] = rgamma(1,(n+a2)/2,(sg+b2)/2)
  }

  return(list(phi1 = phi1, phi2 = phi2, alpha = alph, beta = bet, gamma = gam))
}