Dear R-sig-epi members,
I have a small question concerning a test statistic I would like to construct.
In our study we want to compare the performance of three chromogenic media
(BMR, BR and OX) for detection of Methicillin Resistant Staphylococcus aureus.
Presence of MRSA is checked by PCR.
A typical dataset looks like this:
typ reality plate count
1 bmr yes yes 39
2 bmr yes no 6
3 bmr no yes 2
4 bmr no no 56
5 ox yes yes 40
6 ox yes no 6
7 ox no yes 14
8 ox no no 43
9 br yes yes 43
10 br yes no 5
11 br no yes 14
12 br no no 41
where typ = de type of chromogenic medium
reality = do the bacteria represent MRSA as verified by PCR (yes or no)
plate = are the bacteria recognized as MRSA as observed on the medium? (yes or
no)
or in array format
, , typ = bmr
plate
reality yes no
yes 39 6
no 2 56
, , typ = br
plate
reality yes no
yes 43 5
no 14 41
, , typ = ox
plate
reality yes no
yes 40 6
no 14 43
Thus, for example
typ reality plate count
1 bmr yes yes 39 --> are true positives
3 bmr no yes 2 --> are false positives
In particular, we want to compare the ratio of true positives over true
positives + false negatives among media. For example for bmr this value,
referred to as "spe" = 39/(39+6) = 0.867
To compare the spe's between the three media, I was thinking of a test
statistic comparable to the Breslow-Day statistic (for testing homogeneity of
odds ratios) where our test statistic equals the sum of [(spe[i] -
spe_hat[i])²/spe_hat[i]]. The spe_hat corresponds to the spe value using the
expected frequencies of n11k (n11k_hat) and n12k (n12k_hat)of the i-th partial
table, having the same marginal totals as the observed data, yet having odds
ratio equal to Mantel-Haenszel estimator (the common odds ratio among type
medium).
My problem/question is how to calculate the wanted spe_hat, or else the
expected n11k_hat and n12k_hat ??
Below you find some of my R-code.
Many thanks in advance.
Dieter Anseeuw
mrsa<-read.csv('mrsa.csv', header=T, sep=";")
mrsa.array<-xtabs(count~reality+plate+typ, data=mrsa)
spe.fun<-function(x){
### calculate observed spe's
spe<-x[1,1,]/(x[1,1,]+x[1,2,])
### get estimate of common odds ratio
theta <- mantelhaen.test(x, correct=FALSE)$estimate
### observed frequencies
n11k <- x[1,1,]
n22k <- x[2,2,]
n12k <- x[1,2,]
n21k <- x[2,1,]
### expected frequencies?? spe_hat???
K<-dim(x)[3]
for (i in 1:K){
spe.statistic<-sum((spe[i]-spe_hat[i])^2/spe_hat[i])
}
dfree<-K-1
pval<-1-pchisq(spe.statistic,dfree)
list(test.statistic=spe.statistic, degrees_freedom=dfree, p=pval)
}
--
Dr. Ir. Dieter Anseeuw
Lecturer & Coordinator Afstandsonderwijs Groenmanagement
Katho Campus Roeselare
Wilgenstraat 32
8800 Roeselare Belgium
Direct phone: +32 51 23 29 68
http://www.katho.be/hivb
http://www.linkedin.com/in/dieteranseeuw
[[alternative HTML version deleted]]
_______________________________________________
R-sig-Epi@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-epi
