truncate-platipus.R

stanfile<-write_stan_file(
"data {
  int<lower=2> K;
  int<lower=0> N;
  int<lower=1> D;
  array[N] int<lower=1,upper=K> ybaseline;
    array[N] int<lower=ybaseline, upper=K> y;
  array[N] row_vector[D] x;
}
parameters {
  vector[D] beta;
  ordered[K-1] c;
}

model {
  vector[K] theta;
  real sumtheta;
  for (n in 1:N) {
    real eta;
    eta = x[n] * beta;
    theta[1] = 1 - inv_logit(eta - c[1]);
    for (k in 2:(K-1))
      theta[k] = inv_logit(eta - c[k-1]) - inv_logit(eta - c[k]);
    theta[K] = inv_logit(eta - c[K-1]);
    for(k in 1:(K-1))
       theta[k] = theta[k] * (k >= ybaseline[n]);

    sumtheta = sum(theta);   

     for(k in 1:K)
       theta[k] = theta[k] / sumtheta;
   y[n] ~ categorical(theta);
  }
}"
)


trt<-rbinom(300,1,.5)
latent<-rlogis(300)
manifest<- as.numeric(cut(c(latent,latent-trt/2),15))
outcome<-manifest[1:300]
baseline<-manifest[1:300]
table(trt,outcome)

stan_data <- list(K=15, N=300,  D=1, ybaseline=rep(1,300),y=outcome,x=as.matrix(trt))

mod <- cmdstan_model(stanfile)


fit_mcmc <- mod$sample(
  data = stan_data,
  seed = 123,
  chains = 2,
  parallel_chains = 2
)

fit_mcmc$summary()
MASS::polr(factor(outcome)~trt)

ybaseline<-baseline
ybaseline[outcome==2]<-2
ybaseline[outcome==4]<-4
ybaseline[outcome==9]<-9


stan_data2 <- list(K=15, N=300,  D=1, y=outcome,x=as.matrix(trt), ybaseline=ybaseline)


fit_mcmc2 <- mod$sample(
  data = stan_data2,
  seed = 123,
  chains = 2,
  parallel_chains = 2
)

fit_mcmc2$summary()
MASS::polr(factor(outcome)~trt)