Invalid parent values in estimating an adjacent category logit model

Help
Jun Xu
2014-01-30
2014-01-31
  • Jun Xu

    Jun Xu - 2014-01-30

    Hi all,

    I am trying to estimate an adjacent category logit model. For a quick review, see https://onlinecourses.science.psu.edu/stat504/node/175. I used rjags package in R to estimate the model. I was successful using the model codes at the very end of this email with one data file. But when I switched to another data file, it didn't work any more.

    For the second data file, I got the error message saying that
    "Error: Error in node pr[74,1]
    Invalid parent values"

    In addition, every time after I rerun the codes with the second data file, the node with the error also changes, but it's always about pr[i,1]; that is, the probability of one for some case. I checked the line in the model state related to pr[i,1], and I don't see any inappropriate priors, such as a normal prior for a precision parameter. But maybe I missed something here. For both sets of data, each has a different ordered response variable, ranging from one to four, and each has a different set of independent variables. None of the data sets have missing values. YOu can click the link below this paragraph to download my r script and the two data sets (in Stata format) to reproduce the results and the error. The one data file that worked is gssCum7212BayesOrdSample200.dta (the dependent variable is health) and the other one that didn't work is issp2005Sample200.dta (its dependent variable is worryjob). I sampled from the original data so that the sample sizes are small enough for a quick rerun of my model and data. I just cannot figure out any reason why the codes only worked with some but not other similar data.

    link to my r script and data:

    https://drive.google.com/folderview?id=0Bw7yz4yKob9XVzdrZm5EeG5ESTA&usp=sharing

    Note that you need to change the directory/file name in the r script to run the model off the two data sets. You don't need to change anything else.

    model {
        for( i in 1 : nData ) {
            y[i] ~ dcat( pr[i,1:nYlevels] )
            pr[i,nYlevels] <- 1/denom[i]
            pr[i,3] <- exp(thresh[3]-mu[i])/denom[i]
            pr[i,2] <- exp(thresh[3]-mu[i])*exp(thresh[2]-mu[i])/denom[i]
            pr[i,1] <- exp(thresh[3]-mu[i])*exp(thresh[2]-mu[i])*exp(thresh[1]-mu[i])/denom[i]
            denom[i] <- exp(thresh[3]-mu[i])*exp(thresh[2]-mu[i])*exp(thresh[1]-mu[i])
                        + exp(thresh[3]-mu[i])*exp(thresh[2]-mu[i])
                        + exp(thresh[3]-mu[i])
            mu[i] <- inprod( b[1:nPredictors] , x[i,1:nPredictors] )
        }
        # bPrec <- pow( nYlevels/4 , -2 ) # max plausible slope is 1SD
        for ( j in 1:nPredictors ) {
            b[j] ~ dnorm(0, 0.001) # modest precision because of normalized x,y values
        }
        # threshPriorPrec <- 1
        for ( k in 1:(nYlevels-1) ) {
            threshPriorMean[k] <- k+0.5
            thresh[k] ~ dnorm( threshPriorMean[k] , 0.001)
        }
    }
    
     
  • Matt Denwood

    Matt Denwood - 2014-01-31

    Hi Jun

    Assuming that this error occurs at compilation rather than run time, did you try setting appropriate initial values for the unobserved parents of pr (i.e. b and thresh)? If not, I'd give a named list with a vector (length nPredictors) of zeros for b plus a vector of (1:(nYlevels-1))+0.5 for thresh to the 'inits' argument to jags.model - I suspect this would help with your problem (assuming that these values are actually sensible).

    Hope that helps,

    Matt

     
  • Jun Xu

    Jun Xu - 2014-01-31

    Matt,

    Thanks for your suggestion. Actually I tried feeding some initial values before. I ran a linear regression of y on x (several independent variables), and then I use OLS regression coefficients as the initial values for those coefficients for x. Here is the full r script, which can be downloaded from the link that I provided above:

    require(foreign)
    dataSource = read.dta("data/issp2005Sample200.dta", convert.factor=F)
    
    # adjacent category logit 
    
    require(rjags)
    #------------------------------------------------------------------------------
    # THE MODEL.
    
    modelstring = "
    model {
        for( i in 1 : nData ) {
            y[i] ~ dcat( pr[i,1:nYlevels] )
            pr[i,nYlevels] <- 1/denom[i]
            pr[i,3] <- exp(thresh[3]-mu[i])/denom[i]
            pr[i,2] <- exp(thresh[3]-mu[i])*exp(thresh[2]-mu[i])/denom[i]
            pr[i,1] <- exp(thresh[3]-mu[i])*exp(thresh[2]-mu[i])*exp(thresh[1]-mu[i])/denom[i]
            denom[i] <- exp(thresh[3]-mu[i])*exp(thresh[2]-mu[i])*exp(thresh[1]-mu[i])
                        + exp(thresh[3]-mu[i])*exp(thresh[2]-mu[i])
                        + exp(thresh[3]-mu[i])
            mu[i] <- inprod( b[1:nPredictors] , x[i,1:nPredictors] )
        }
        # bPrec <- pow( nYlevels/4 , -2 ) 
        for ( j in 1:nPredictors ) {
            b[j] ~ dnorm(0, 0.001) 
        }
        # threshPriorPrec <- 1
        for ( k in 1:(nYlevels-1) ) {
            threshPriorMean[k] <- k+0.5
            thresh[k] ~ dnorm( threshPriorMean[k] , 0.001)
        }
    }
    " # close quote for modelstring
    writeLines(modelstring,con="model.txt")
    
    #------------------------------------------------------------------------------
    # THE DATA.
    
    dataMat = as.matrix(dataSource)
    
    # Rename for use by generic processing later:
    nData = NROW(dataMat)
    x = dataMat[,-1]
    predictorNames = colnames(dataMat)[-1]
    nPredictors = NCOL(x)
    y = as.matrix(dataMat[,1])
    predictedName = colnames(dataMat)[1]
    nYlevels = max(y)
    
    lmInfo = lm( y ~ x ) # R function returns MLE
    # bInit = lmInfo$coef[-1]
    bInit = rep(0, nPredictors)
    
    dataList = list(
      x = x ,
      y = as.vector( y ) , # BUGS does not treat 1-column mat as vector
      nPredictors = nPredictors ,
      nData = nData ,
      nYlevels = nYlevels
    )
    
    #------------------------------------------------------------------------------
    # INTIALIZE THE CHAINS.
    
    initsList = list(
      b = bInit , # from lm(y~x), above
      thresh = 1:(nYlevels-1)+.5
    )
    
    #------------------------------------------------------------------------------
    # RUN THE CHAINS
    
    parameters = c( "b" , "thresh" )  
    adaptSteps = 500              # Number of steps to "tune" the samplers.
    burnInSteps = 500            # Number of steps to "burn-in" the samplers.
    nChains = 3                   # Number of chains to run.
    numSavedSteps=50000           # Total number of steps in chains to save.
    thinSteps=1                   # Number of steps to "thin" (1=keep every step).
    nPerChain = ceiling( ( numSavedSteps * thinSteps ) / nChains ) # Steps per chain.
    # Create, initialize, and adapt the model:
    jagsModel = jags.model( "model.txt" , data=dataList , inits=initsList , 
                            n.chains=nChains , n.adapt=adaptSteps )
    # Burn-in:
    cat( "Burning in the MCMC chain...\n" )
    update( jagsModel , n.iter=burnInSteps )
    # The saved MCMC chain:
    cat( "Sampling final MCMC chain...\n" )
    codaSamples = coda.samples( jagsModel , variable.names=parameters , 
                                n.iter=nPerChain , thin=thinSteps )
    
    mcmcChain = as.matrix( codaSamples )
    gelman.diag(codaSamples)
    
     
  • Matt Denwood

    Matt Denwood - 2014-01-31

    It seems your model is incredibly sensitive to the value of b[] - I can get it to run by constraining that parameter to within -0.1,0.1 (change b[j] ~ dnorm(0, 0.001) to b[j] ~ dnorm(0, 0.001)T(-0.1,0.1) in your model). Any wider range consistently gives the error you noted, presumably caused by sampling values of b outside this range which consequently gives extremely implausible values for pr (the form of the equation for pr makes me think it may be numerically quite unstable, which is one possible explanation). If it should really be possible for b to exist outside this range (particularly if the posterior for b goes anywhere near the limits of the constrained parameter space) then you should look at re-formulating your model - others may be able to help you with this, but I'm afraid I'm not particularly familiar with this class of model.

    Hope that helps,

    Matt

     
    Last edit: Matt Denwood 2014-01-31
  • Jun Xu

    Jun Xu - 2014-01-31

    By incorporating your suggested change, my codes now worked. With that in my mind and the fact the my initial codes did work for another data set, I start to poke around the second data file that didn't work with my codes. I noticed that the values for age and especially agesq are quite large. So I rescaled these two variables, and it turns out the my codes would've worked had the x's been properly value. So thanks a lot. By the way, what is the T function does? Does it confine the final sampling to be between certain -0.1 and 0.1? Is this like truncated normal (sampling only from -0.1 to 0.1?

     
    Last edit: Jun Xu 2014-01-31
  • Matt Denwood

    Matt Denwood - 2014-01-31

    Glad that helped get you onto the right track - re-scaling (and also possibly re-centring?) of the predictors sounds like a good idea :)

    Censoring/truncation of variables is detailed (far better than I could) in section 7.0.6 of the JAGS user manual. If you are keen to think a little bit deeper about it, then have a read of this post on Martyn's blog: http://martynplummer.wordpress.com/2010/03/23/unknown-unknowns/

     
  • Jun Xu

    Jun Xu - 2014-01-31

    This is wonderful. I am relatively new to Jags, and only replicated some Jags codes from sources on related likelihood-based models. It's important to dig into these nuances. Thanks a lot!

     

Get latest updates about Open Source Projects, Conferences and News.

Sign up for the SourceForge newsletter:





No, thanks