Add effects to difference between cutpoints in multilevel ordered logistic model

Modeling
1

I am building an ordered logistic model with three categories. (Trees in phenological states determined by forcing temperatures). I would like to model the difference between cutpoints with effects from different clusters in my data.

I’ve created the transformed parameter

cdiff = cutpoints[2]-cutpoints[1];

and put a prior on that parameter.

 alpha ~ lognormal(1,1);
 cdiff ~ gamma(alpha, 1); // difference between first and second cutpoint

I can include the difference between cutpoints, but I can’t figure out how to include the clusters in the difference.

Some of the clusters are crossed and others are nested. (It’s a messy common garden design). Clones are represented by individual trees, which are male or female. Clones have a provenance. Provenances and clones are both represented at multiple sites, but an individual tree can obviously only occur at one site).

What I’d like to do is something like

alpha = \alpha + \alpha_{site} + \alpha_{provenance} +\alpha_{clone\ in\ provenance} + \alpha_{year}

I am new to stan and can’t figure out how to add effects to the difference between cutpoints.

Here is a working version for simulated data with no groups.

Simulate data in R and run model

cuts <- c(15,20)
beta <- 0.05
forcing <- runif(300, 175, 500) #                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         
eta <- beta * forcing
state <- rordlogit(length(eta), phi = eta, a = cuts)
sim <- data.frame(forcing,state)

K = length(unique(sim$state))
N = nrow(sim)
stan_rdump(c("N", "K", "forcing", "state"), "simulated_nogroups.Rdump")
simplesimdat <- read_rdump("simulated_nogroups.Rdump")

simfit <- stan("cutpoints_difference.stan",
               chains=5, cores=5, iter=2000,
               data=simplesimdat)

Stan code

data{
    int N; //number of observations
    int K; //number of possible states
    int state[N];
    vector[N] forcing; //predictor
}

parameters{
    positive_ordered[K-1] cutpoints;
    real<lower=0,upper=1> beta; // eta "slope"
    real<lower=0> alpha;
}

transformed parameters {
    //declare params
    vector[K-1] h50; // transition points on forcing scale
    real<lower=0> cdiff; //difference between cutpoints
    //define params
    for ( i in 1:2 ) {
        h50[i] = cutpoints[i]/beta;
    }
    cdiff = cutpoints[2]-cutpoints[1];
}

model{
    //declarations
    vector[N] eta;
    //priors
        //for eta
    beta ~ beta( 0.5 , 5 );
        //for cutpoints
    cutpoints[1] ~ exponential(0.3); // don't let first cutpoint get too big for its britches
    alpha ~ lognormal(1,1);
    cdiff ~ gamma(alpha, 1); // difference between first and second cutpoint
    //model
    for ( i in 1:N ) {
        eta[i] = beta * forcing[i];
    }
        state ~ ordered_logistic( eta , cutpoints );
}

I’m not actually interested in the cutpoints, I’m interested in the point at which half of the trees have transitioned from one state to another (h50), defined in the transformed parameters as cutpoints/beta.

I am trying to use the difference between because I believe that some effects should shift the h50 in different directions. If I include the effects in eta, then h50 can only be shifted in one direction by an effect for both transitions.

I am also open to adding effects to both cutpoints. I tried to include effects on each cutpoint individually based on this question. However I couldn’t figure out how to include multiple clusters and was only able to include one cluster.

Is it possible to add effects to the difference between cutpoints? How can I do this?

2

I don’t understanding your problem fully, but some suggestions that may or may not be useful:

  • AFAIK people have cutpoints[1] fixed to 0 or 1 in most ordered logistic models as it is not necessarily identified (but not sure this applies to your case)
  • Once you have cutpoints[1] you can compute cdiff from your predictors. Since cdiff needs to be positive, a natural way would be to exponentiate the predictor cdiff = exp(cdiff_predictor)

Hope that helps.

1 Like