Calculating Kappa Coefficient of Agreement

This example illustrates hows JAGS handles a more complex hierarchical Bayesian model and compares the output from JAGS with the output from WinBUGS

The model was taken from the WinBUGS course developed by Michael Lee and Eric-Jan Wagenmakers at http://www.ejwagenmakers.com/BayesCourse/BayesBook.html

Suppose we have data in the form of agreement counts between two methods, the objective method and the surrogate method. The data is represented by a vector d = [ a b c d ], where

The goal is to calculate Kappa, a standard measure of agreement

Contents

Model definition

We use a hierarchical Bayesian model to estimate Kappa:

# Kappa Coefficient of Agreement
model {
  # Underlying Rates
  # Rate objective method decides 'one'
  alpha ~ dbeta(1,1)
  # Rate surrogate method decides 'one' when objective method decides 'one'
  beta ~ dbeta(1,1)
  # Rate surrogate method decides 'zero' when objective method decides 'zero'
  gamma ~ dbeta(1,1)
  # Probabilities For Each Count
  pi[1] <- alpha*beta
  pi[2] <- alpha*(1-beta)
  pi[3] <- (1-alpha)*(1-gamma)
  pi[4] <- (1-alpha)*gamma
  # Count Data
  d[1:4] ~ dmulti(pi[],n)
  # Derived Measures
  # Rate surrogate method agrees with the objective method
  xi <- alpha*beta+(1-alpha)*gamma
  # Rate of chance agreement
  psi <- (pi[1]+pi[2])*(pi[1]+pi[3])+(pi[2]+pi[4])*(pi[3]+pi[4])
  # Chance corrected agreement
  kappa <- (xi-psi)/(1-psi)
}

clear;
clc;

Defining Observed Data

d=[14 4 5 210];
%d=[20 7 103 417];
%d=[157 0 13 0];

% Derived Variables
n=sum(d);

% Create single Matlab structure with values for all the observed JAGS nodes
datastruct = struct('d',d,'n',n);

% MCMC parameters for JAGS
nchains=7; % How Many Chains?
nburnin=500; % How Many Burn-in Samples?
nsamples=1000; % How Many Recorded Samples?

% Initialize values all latent variables in all chains
for i=1:nchains
    S.alpha = 0.5; % An Initial Value
    S.beta = 0.5; % An Initial Value
    S.gamma = 0.5; % An Initial Value
    init0(i) = S;
end

Use JAGS to Sample

tic
doparallel = 0;
fprintf( 'Running JAGS\n' );
[samples, stats ] = matjags( ...
        datastruct, ...
        fullfile(pwd, 'Kappa_1.txt'), ...
        init0, ...
        'doparallel' , doparallel, ...
        'nchains', nchains,...
        'nburnin', nburnin,...
        'nsamples', nsamples, ...
        'thin', 1, ...
        'monitorparams', {'kappa','xi','psi','alpha','beta','gamma','pi'}, ...
        'savejagsoutput' , 1 , ...
        'verbosity' , 1 , ...
        'cleanup' , 0  );
Running JAGS
Running chain 1 (serial execution)
Running chain 2 (serial execution)
Running chain 3 (serial execution)
Running chain 4 (serial execution)
Running chain 5 (serial execution)
Running chain 6 (serial execution)
Running chain 7 (serial execution)
toc
Elapsed time is 1.732531 seconds.

Get the posterior means from JAGS:

stats.mean
ans = 

    kappa: 0.6975
       xi: 0.9544
      psi: 0.8474
    alpha: 0.0811
     beta: 0.7488
    gamma: 0.9726
       pi: [0.0607 0.0204 0.0252 0.8937]

Use WinBUGS to Sample

tic
fprintf( 'Running WinBUGS\n' );
[samples2, stats2] = matbugs(datastruct, ...
    fullfile(pwd, 'Kappa_1.txt'), ...
    'init', init0, ...
    'nChains', nchains, ...
    'view', 0, 'nburnin', nburnin, 'nsamples', nsamples, ...
    'thin', 1, 'DICstatus', 0, 'refreshrate',100, ...
    'monitorParams', {'kappa','xi','psi','alpha','beta','gamma','pi'}, ...
    'Bugdir', 'C:/WinBUGS14');
toc
Running WinBUGS
Elapsed time is 4.273773 seconds.

Get the posterior means from WinBUGS:

stats2.mean
ans = 

    alpha: 0.0810
     beta: 0.7505
    gamma: 0.9726
    kappa: 0.6986
       pi: [0.0608 0.0202 0.0252 0.8938]
      psi: 0.8474
       xi: 0.9546