### Simple simulation – Am I being naive?

Saw this sometime last week: http://bengoldacre.posterous.com/i-think-this-guardian-story-is-a-bit-wrong

I will not rehash the main points, but the idea that popped up in my mind was, the patients will not necessarily all have the same mortality rate so what happens if we have a population of patients with different mortality rates?

I thought of possible ways to simulate this, but did not sit down to try it till today. Took me a while to figure out the relevant MATLAB commands and to structure it. Thanks goes to Leonard who gave me some ideas on how to write the code.

Simple assumptions of population of patients:

- 10% have 6% mortality rate (difficult case)
- 25% have 3% mortality rate (moderate case)
- 65% have 1% mortality rate (simple case)

And then each hospital has a certain number of patients, which is fixed. How do the hospitals do?

For size 5:

Value Count Percent

0 2653 26.53%

1 5081 50.81%

2 1946 19.46%

3 301 3.01%

4 19 0.19%

For size 10:

Value Count Percent

0 541 5.41%

1 2118 21.18%

2 3224 32.24%

3 2488 24.88%

4 1188 11.88%

5 370 3.70%

6 67 0.67%

7 4 0.04%

I did 10000 realizations on both.

Would it be fair then to say that a hospital with low number of patients can easily get a high mortality rate by just being unlucky?

MATLAB code for those interested below. I may have missed out something as I finished this in less than an hour!

-------------------------------------------------------------------------------

function AAAsim(nRlz,sz)

% Simulation of AAA patients.

% Assume 2% mortality in general

% 3 Classes of Patients:

% 10% chance of patient with 6% mortality

% 25% chance of patient with 3% mortality

% 65% chance of patient with 1% mortality

%

% Assume each hospital has sz patients

n = nRlz * sz; %population size

index = randperm(n); % Randomly order patients

for i=1:nRlz

a = sz*(nRlz - 1) + 1;

b = sz*(nRlz);

sample = index(a:b);

sam_death = 0;

for j=1:sz

pt = sample(j);

if(pt <=(0.1*n));

prob = 0.6;

elseif(pt > (0.1*n) && pt <=(0.35*n));

prob = 0.3;

else(pt > (0.35*n));

prob = 0.1;

end

prob;

ptdoa = binornd(1,prob); % Returns 0 if alive, 1 if dead

sam_death = sam_death + ptdoa;

end

deaths(i)= sam_death;

fprintf('Iteration number %d completed.\n',i);

end

tabulate(deaths)

hist(deaths,0:sz)

## No comments

## Add a comment

You are not allowed to comment on this entry as it has restricted commenting permissions.