Kappa Cohen Estimator

Goodevening everyone.I am trying to estimate Kappa Cohen coefficient with Rstan but something is wrong.Can you please help me.I am new in Stan.


# Kappa Cohen 

KappaCohen <- "
data {
  int nA;
  int nB;
  int nC;
  int nD;
}

parameters {
  real alpha;
  real beta;
  real gamma;
  real<lower=0, upper=1> pA;
  real<lower=0, upper=1> pB;
  real<lower=0, upper=1> pC;
  real<lower=0, upper=1> pD;
  
}

model {
  alpha ~ beta(1, 1);
  beta  ~ beta(1, 1);
  gamma ~ beta(1, 1);
}

generated quantities {
  pA = alpha*beta ;
  pB = alpha*(1-beta) ;
  pC =(1-alpha)*(1-gamma) ; 
  pD = (1-alpha)*gamma ; 
  real<lower=0> xi;
  real<lower=0> psi;
  real<lower=0> kappa;
  xi = alpha * beta + (1-alpha)*gamma ; 
  psi = (pA + pB) *(pA + pC) + (pB+pD)*(pC+pD);
  kappa = (xi - psi)/(1-psi);
}
"
data_list <- list(nA = 14, nB = 4, nC = 5, nD = 210)

stan_samplesKappa <- stan(model_code = KappaCohen, data = data_list)

You really should show the error messages in your output instead of just saying that “something” is wrong.

Anyway, at least one thing that I notice is that pA, pB, pC, and pD are not used in the model block, but are listed as parameters. They should be declared in the “generated quantities” block instead. Also, variable declarations have to be the first statements in a block.

Got it


# Kappa Cohen 
KappaCohen <- "
// Kappa Coefficient of Agreement
data { 
  int<lower=0> y[4];   
}
parameters {
  // Underlying Rates
  // Rate first rater  decides 'one'
  real<lower=0,upper=1> alpha;
  // Rate Ssecond rater decides 'one' when the first rater  decides 'one'
  real<lower=0,upper=1> beta;
  // Rate second rater  decides 'zero' when the first rater decides 'zero'
  real<lower=0,upper=1> gamma;
} 
transformed parameters {
  simplex[4] pi;
  real xi;
  real psi;
  real kappa;
  // Probabilities For Each Count
  pi[1] = alpha * beta;
  pi[2] = alpha * (1 - beta);
  pi[3] = (1 - alpha) * (1 - gamma);
  pi[4] = (1 - alpha) * gamma;
    
  
  
  xi = alpha * beta + (1 - alpha) * gamma ;
  psi = (pi[1] + pi[2]) * (pi[1] + pi[3]) + (pi[2] + pi[4]) * (pi[3] + pi[4]);  
  kappa = (xi - psi) / (1 - psi);   
}
model {
  alpha ~ beta(1, 1);  
  beta ~ beta(1, 1);  
  gamma ~ beta(1, 1);  
  y ~ multinomial(pi);   // Count Data 
}"



y <- c(14, 4, 5, 210)