Dealing with Autocorrelation

2013-02-25
2013-02-27
  • mark pahuta

    mark pahuta - 2013-02-25

    Hi,

    I am using dmstate to fit a multistate model. The data are collected from several different studies, so I have incorporated a random effects term. I am also testing a frailty model.

    My baseline transition intensity is a[,]~dgamma(1E-4, 1E-4) T (1E-4,). My time-scale is days, and median time is around 2 years. I've also experimented with a uniform prior on log(2)/a~dunif.

    I am having extreme autocorrelation for a, beta is OK. The model takes a very long time to run, so I was hoping to try something other than thinning. The model adapts successfully. I hve autocorrelation with the gamma and uniform prior.

    My code is below. Does anyone have suggestions on other things I could do to reduce autocorrelation?

    model{
     for(pat in 1:Npat){
      So[pat,1]~dcat(impstates[1,,Pat.CovarCombo.Study[pat,3]])
      for (obs in 2:nobs[pat]){
          censtyp[pat,obs] ~ dinterval(So[pat,obs], brkpt[pat,obs,])  
          So[pat,obs] ~ dmstate (So[pat,obs-1], tm_int[pat,obs], Q[,,pat,piece[pat,obs]])
      }
      }
    
    for(study in 1:3){
    impstates[1,1,study]~dunif(0.025,0.975)
    impstates[1,2,study]<-1-impstates[1,1,study]
    }
    
    #Frailty and Qmatrix For Each Participant
    for (pat in 1:Npat){
    frail[pat]~dnorm(0,sigma.frail)
    
    for (pieceid in 1:maxpiecepat[pat]){
    
    Q[1,2,pat,pieceid] <- a[1,pieceid]*exp(alpha[1,Pat.CovarCombo.Study[pat,3]]+theta[Pat.CovarCombo.Study[pat,2],1]+frail[pat])
    Q[1,3,pat,pieceid] <- a[2,pieceid]*exp(alpha[2,Pat.CovarCombo.Study[pat,3]]+theta[Pat.CovarCombo.Study[pat,2],2]+frail[pat])
    Q[2,3,pat,pieceid] <- a[3,pieceid]*exp(alpha[3,Pat.CovarCombo.Study[pat,3]]+theta[Pat.CovarCombo.Study[pat,2],3]+frail[pat])
    
    Q[1,1,pat,pieceid]<- -1 * (Q[1,2,pat,pieceid] + Q[1,3,pat,pieceid])  
    Q[2,1,pat,pieceid]<-0
    Q[2,2,pat,pieceid]<- -1 * Q[2,3,pat,pieceid]
    Q[3,1,pat,pieceid]<-0
    Q[3,2,pat,pieceid]<-0
    Q[3,3,pat,pieceid]<-0
    }}
    
    #Covar x Coef Combos
    for(num_cov in 1:num_covar_combo){
    for(trans in 1:3){
    theta[num_cov,trans] <- inprod(UniqueCovarCombo[num_cov,], beta[trans,])
    }}
    
    #Coefficient For Each Transition
    for (trans in 1:3){
    for (covar in 1:num_covars){
    beta[trans,covar]~dnorm(0, sigma.beta)
    }}
    
    #Random Effect For Each Study
    for (study in 1:3){
    alpha[1:3,study] ~ dmnorm(mu.alpha[1:3], Omega.alpha[1:3,1:3,1])
    }
    
    #Prior For Each Transition Intensity
    for (trans in 1:3){
    for (pieceid in 1:max_mod_piece){
    a[trans,pieceid]~dgamma(1E-4, 1E-4) T (1E-4,)
    }}
    
    Omega.alpha[1:3,1:3,1] ~ dwish(Lambda, nu )
    #Scale Prior For Each Transition Intensity
    Lambda[1,1]<-Order.R.mat
    Lambda[1,2]<-0
    Lambda[1,3]<-0
    Lambda[2,1]<-0
    Lambda[2,2]<-Order.R.mat
    Lambda[2,3]<-0
    Lambda[3,1]<-0
    Lambda[3,2]<-0
    Lambda[3,3]<-Order.R.mat
    }
    
     
    Last edit: mark pahuta 2013-02-25
  • Luc Coffeng

    Luc Coffeng - 2013-02-27

    Have you tried using redundant parameters? E.g., specify a parameter in the model as the sum or product of two 'raw' parameters? I'm not sure how this would exactly apply to your model, but perhaps the following example of a simple random effects may provide inspiration. In this example, the parameters xi1 and xi2 are redundant.

    model {
        for (i in 1:N) {
            # Likelihood of type1 binomially distributed data
            k1[i] ~ dbin(p1[i],n1[i])
            logit(p1[i]) <- B0[1] + e.vill[i,1]
    
            # Likelihood of type2 binomially distributed data
            k2[i] ~ dbin(p2[i],n2[i])
            logit(p2[i]) <- B0[2] + e.vill[i,2]
    
            # Correlated random effects
            e.vill[i,1] <- e.vill.raw[i,1] * xi1
            e.vill[i,2] <- e.vill.raw[i,2] * xi2
            e.vill.raw[i,1:2] ~ dmnorm(Mu,Sigma2.inv.raw) # Mu = c(0,0)
        }
        # Scaled inverse Wishart prior for random village effects
        Sigma2.inv.raw ~ dwish(R,3) # R is a 2 by 2 identity matrix
        xi1 ~ dunif(0,100)
        xi2 ~ dunif(0,100)
        Sigma2.raw <- inverse(Sigma2.inv.raw)
        sigma1 <- pow(Sigma2.raw[1,1],0.5) * xi1
        sigma2 <- pow(Sigma2.raw[2,2],0.5) * xi2
        rho12 <- Sigma2.raw[1,2]/sqrt(Sigma2.raw[1,1]*Sigma2.raw[2,2])
    }
    

    There are more ways to perform redundant parameterization, and you have to try which is best for your model. Andrew Gelman's book Data Analysis Using Regression and Multilevel/Hierarchical Models provides helpful examples. Hope this helps!

    Luc

     
    Last edit: Luc Coffeng 2013-02-27

Log in to post a comment.