Return to Question

I want to model the index called Delta P (e.g., p.144 of this paper), which is basically a difference between two proportions (i.e., $\frac{n_1}{N_1}$ - $\frac{n_2}{N_2}$), as a function of a predictor. The input data should be the four count variables from which to calculate Delta P (i.e., $n_1$, $N_1$, $n_2$, $N_2$) and predictor values.

[Note on cross-posting: This question has now been posted on the Stan Forums as well.]

I want to model the index called Delta P (e.g., p.144 of this paper), which is basically a difference between two proportions (i.e., $\frac{n_1}{N_1}$ - $\frac{n_2}{N_2}$), as a function of a predictor. The input data should be the four count variables from which to calculate Delta P (i.e., $n_1$, $N_1$, $n_2$, $N_2$) and predictor values.

Below is my attempt to do it in Stan. When I run the code, I get a message about Jacobian adjustments since the left-hand side of a sampling statement is deltaP, which is calculated by subtracting two parameters (theta1 and theta2, which are the estimated values of $\frac{n_1}{N_1}$ and $\frac{n_2}{N_2}$, respectively).

data { 
  int<lower=0> N; // total number of observations
  int<lower=1> denom1[N]; // denominator of the first proportion
  int<lower=1> denom2[N]; // denominator of the second proportion
  int<lower=0> nom1[N]; // nominator of the first proportion
  int<lower=0> nom2[N]; // nominator of the second proportion
  real x[N]; // predictor variable
} 

parameters {
  real<lower=0, upper=1> theta1[N]; // the first proportion
  real<lower=0, upper=1> theta2[N]; // the second proportion
  real alpha; // intercept
  real beta; // slope parameter for x
  real<lower=0> sigma; // SD of the error term
} 

transformed parameters {
  real<lower=-1, upper=1> deltaP[N]; // Delta P
  for (i in 1:N) {
    deltaP[i] = theta1[i] - theta2[i];
  }
}

model {
  // priors
  theta1 ~ beta(1, 1);
  theta2 ~ beta(1, 1);
  alpha ~ normal(0, 2);
  beta ~ normal(0, 2);
  sigma ~ normal(0, 1) T[0, ];
  
  for (i in 1:N) {
    // estimating thetas based on denoms and noms
    nom1[i] ~ binomial(denom1[i], theta1[i]); 
    nom2[i] ~ binomial(denom2[i], theta2[i]);
    // deltaP is sampled from the truncated normal distribution whose mean is alpha + beta * x and the SD is sigma
    deltaP[i] ~ normal(alpha + beta * x[i], sigma) T[-1, 1];
  }
}

library("rstan")

### Generate fake data
set.seed(100)
# sample size
N <- 100
# True parameter values
alpha <- -0.2
beta <- 0.5
sigma <- 0.1

# predictor values (x) and Delta P values
while (TRUE) {
  x <- runif(N, -1, 1)
  deltaP <- alpha + beta * x + rnorm(N, sd = sigma)
  if (all(deltaP <= 1) & all(deltaP >= -1)) break
}
# theta values
theta1 <- theta2 <- numeric(N)
for (i in 1:N) {
  if (deltaP[i] > 0) {
    theta1[i] <- runif(1, deltaP[i], 1)
    theta2[i] <- theta1[i] - deltaP[i]
  } else {
    theta2[i] <- runif(1, abs(deltaP[i]), 1)
    theta1[i] <- theta2[i] + deltaP[i]
  }
}

# denoms and noms
denom1 <- sample(N, replace = TRUE)
denom2 <- sample(N, replace = TRUE)
nom1 <- rbinom(N, denom1, theta1)
nom2 <- rbinom(N, denom2, theta2)

### fit the model
fit <- stan(file = 'xxx.stan', 
            data = list(
              N = N,
              denom1 = denom1,
              denom2 = denom2,
              nom1 = nom1,
              nom2 = nom2,
              x = x
            ))

DIAGNOSTIC(S) FROM PARSER:
Info:
Left-hand side of sampling statement (~) may contain a non-linear transform of a parameter or local variable.
If it does, you need to include a target += statement with the log absolute determinant of the Jacobian of the transform.
Left-hand-side of sampling statement:
    deltaP[i] ~ normal(...)

Below is my attempt to do it in Stan. When I run the code, I get a message about Jacobian adjustments since the left-hand side of a sampling statement is deltaP, which is calculated by subtracting one parameter from another (theta1 - theta2, where theta1 is the estimated value of $\frac{n_1}{N_1}$ and theta2 is that of $\frac{n_2}{N_2}$).

data { 
  int<lower=0> N; // total number of observations
  int<lower=1> denom1[N]; // denominator of the first proportion
  int<lower=1> denom2[N]; // denominator of the second proportion
  int<lower=0> nom1[N]; // nominator of the first proportion
  int<lower=0> nom2[N]; // nominator of the second proportion
  real x[N]; // predictor variable
} 

parameters {
  real<lower=0, upper=1> theta1[N]; // the first proportion
  real<lower=0, upper=1> theta2[N]; // the second proportion
  real alpha; // intercept
  real beta; // slope parameter for x
  real<lower=0> sigma; // SD of the error term
} 

transformed parameters {
  real<lower=-1, upper=1> deltaP[N]; // Delta P
  for (i in 1:N) {
    deltaP[i] = theta1[i] - theta2[i];
  }
}

model {
  // priors
  theta1 ~ beta(1, 1);
  theta2 ~ beta(1, 1);
  alpha ~ normal(0, 2);
  beta ~ normal(0, 2);
  sigma ~ normal(0, 1) T[0, ];
  
  for (i in 1:N) {
    // estimating thetas based on denoms and noms
    nom1[i] ~ binomial(denom1[i], theta1[i]); 
    nom2[i] ~ binomial(denom2[i], theta2[i]);
    // deltaP is sampled from the truncated normal distribution whose mean is alpha + beta * x and the SD is sigma
    deltaP[i] ~ normal(alpha + beta * x[i], sigma) T[-1, 1];
  }
}

library("rstan")

### Generate fake data
set.seed(100)
# sample size
N <- 100
# True parameter values
alpha <- -0.2
beta <- 0.5
sigma <- 0.1

# predictor values (x) and Delta P values
while (TRUE) {
  x <- runif(N, -1, 1)
  deltaP <- alpha + beta * x + rnorm(N, sd = sigma)
  if (all(deltaP <= 1) & all(deltaP >= -1)) break
}
# theta values
theta1 <- theta2 <- numeric(N)
for (i in 1:N) {
  if (deltaP[i] > 0) {
    theta1[i] <- runif(1, deltaP[i], 1)
    theta2[i] <- theta1[i] - deltaP[i]
  } else {
    theta2[i] <- runif(1, abs(deltaP[i]), 1)
    theta1[i] <- theta2[i] + deltaP[i]
  }
}

# denoms and noms
denom1 <- sample(N, replace = TRUE)
denom2 <- sample(N, replace = TRUE)
nom1 <- rbinom(N, denom1, theta1)
nom2 <- rbinom(N, denom2, theta2)

### fit the model
fit <- stan(file = 'xxx.stan', 
            data = list(
              N = N,
              denom1 = denom1,
              denom2 = denom2,
              nom1 = nom1,
              nom2 = nom2,
              x = x
            ))

DIAGNOSTIC(S) FROM PARSER:
Info:
Left-hand side of sampling statement (~) may contain a non-linear transform of a parameter or local variable.
If it does, you need to include a target += statement with the log absolute determinant of the Jacobian of the transform.
Left-hand-side of sampling statement:
    deltaP[i] ~ normal(...)