# 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)
``````