[R] OT: A test with dependent samples.

Rolf Turner r.turner at auckland.ac.nz
Wed Feb 11 19:51:22 CET 2009

Thank you ***VERY*** much.  That clarifies a lot of things which were  
knocking around in my mind but about which I could not get my  
thoughts properly

Thanks again.


		Rolf Turner

On 11/02/2009, at 11:27 PM,  
Meyners,Michael,LAUSANNE,AppliedMathematics wrote:

> Rolf,
> as you explicitly asked for a comment on your proposal: It is  
> generally
> equivalent to McNemar's test and maybe even more appropriate  
> because of
> the asymptotics involved in the chi-squared distribution, which might
> not be too good with n=12. In some more detail:
> McNemar's test basically considers the difference between the
> off-diagonal elements, normalized by their sum. You do the same,
> ignoring all (0,0) and (1,1) pairs (the latter do not occur in your
> example anyway). binom.test(12,12,alternative="greater") gives you the
> one-sided p-value.
> If, however, you do not use the exact binomial distribution (with  
> n=12)
> but the normal approximation, you find E(X)=6 and VAR(X)=3 with X the
> number of sequences (0,1), and you observed x=12. This gives you a
> z-score of (12-6)/sqrt(3) and a corresponding one-sided p-value:
>> pnorm(6/sqrt(3), lower.tail=F)
> [1] 0.0002660028
> Further, McNemar's test WITHOUT continuity correction gives
>> mcnemar.test(matrix(c(61,0,12,0),2,2), correct=F)
>         McNemar's Chi-squared test
> data:  matrix(c(61, 0, 12, 0), 2, 2)
> McNemar's chi-squared = 12, df = 1, p-value = 0.000532
> It comes as no surprise that the reported (two-sided!) p-value is
> exactly twice the (one-sided!) p-value from the binomial test with
> normal approximation:
>> pnorm(6/sqrt(3), lower.tail=F)*2
> [1] 0.0005320055
> I do not want to stress the pros and cons of continuity  
> corrections, but
> if neglected, you see that you get the same results (except that  
> McNemar
> is generally two-sided), due to the relation of normal and chi-squared
> distribution.
> If you use the binomial test, you can forget about asymptotics and
> continuity correction, that's why I indeed consider you approach
> superior to the chi-squared approximation used in McNemar's test. But
> you might fail to find a reference for the exact approach...
> You briefly asked in a later mail testing for p=0. Indeed, _any_
> incident will disprove this hypothesis, and the p value reported by
>> binom.test(1,73,p=0)
>         Exact binomial test
> data:  1 and 73
> number of successes = 1, number of trials = 73, p-value < 2.2e-16
> alternative hypothesis: true probability of success is not equal to 0
> 95 percent confidence interval:
>  0.0003467592 0.0739763232
> sample estimates:
> probability of success
>             0.01369863
> is not wrong (as R just tells us it is BELOW a certain value), but  
> could
> be refined by saying "p-value = 0". BTW, I am not sure what
> "p-value=TRUE" tells me, as derived from binom.test(0,73,p=0). I
> personally don't care about either, as to test H0: p=0, I would not  
> use
> any software but rather stress some basic statistical theory.
> It remains the questions whether McNemar (or your proposal) is  
> generally
> appropriate here. Like others, I have doubts, except maybe if
> observations were made on two subsequent days, on the second of which
> the drug was administered. How could you otherwise be sure that the
> increase in incidences it not due to the progression of cancer, say,
> which would be difficult to rule out. Also, from my experience and in
> contrast to you point of view, I'd rather assume it likely that a  
> vet is
> much more reluctant to apply the new treatment to an already vomitting
> cat, your "baseline" with 0 out of 73 at least seems suspicious, and
> until proven wrong, I'd take the assumption that this is not by chance
> only.
> I understand that this is an observational study with no
> control/randomization, but concluding a side-effect from a  
> statistically
> significant change in incidence rates from this study seems  
> questionable
> to me. What I'd propose (and used in comparable situations) is to
> confine yourself on the confidence interval and let the investigator
> decide / discuss whether the lower bound of this is higher than  
> what she
> would expect under control or alternative treatment. She needs to have
> some experience, if not historical data, and using a test or just this
> confidence interval, you will only get a hint on a possible side- 
> effect,
> no formal proof (which is always difficult from observational  
> studies).
> Discussing the CI would seem fairer and even stronger to me then. In
> other words: Mathematically, you'll get an appropriate test, while
> conceptually, I'm far from being convinced.
> Hope this makes sense,
> Michael
> -----Original Message-----
> From: r-help-bounces at r-project.org [mailto:r-help-bounces at r- 
> project.org]
> On Behalf Of Rolf Turner
> Sent: Dienstag, 10. Februar 2009 22:33
> To: R-help Forum
> Subject: [R] OT: A test with dependent samples.
> I am appealing to the general collective wisdom of this list in  
> respect
> of a statistics (rather than R) question.  This question comes to me
> from a friend who is a veterinary oncologist.  In a study that she is
> writing up there were 73 cats who were treated with a drug called
> piroxicam.  None of the cats were observed to be subject to vomiting
> prior to treatment; 12 of the cats were subject to vomiting after
> treatment commenced.  She wants to be able to say that the  
> treatment had
> a ``significant''
> impact with respect to this unwanted side-effect.
> Initially she did a chi-squared test.  (Presumably on the matrix
> matrix(c(73,0,61,12),2,2) --- she didn't give details and I didn't
> pursue
> this.) I pointed out to her that because of the dependence --- same 73
> cats pre- and post- treatment --- the chi-squared test is  
> inappropriate.
> So what *is* appropriate?  There is a dependence structure of some  
> sort,
> but it seems to me to be impossible to estimate.
> After mulling it over for a long while (I'm slow!) I decided that a
> non-parametric approach, along the following lines, makes sense:
> We have 73 independent pairs of outcomes (a,b) where a or b is 0 if  
> the
> cat didn't barf, and is 1 if it did barf.
> We actually observe 61 (0,0) pairs and 12 (0,1) pairs.
> If there is no effect from the piroxicam, then (0,1) and (1,0) are
> equally likely.  So given that the outcome is in {(0,1),(1,0)} the
> probability of each is 1/2.
> Thus we have a sequence of 12 (0,1)-s where (under the null  
> hypothesis)
> the probability of each entry is 1/2.  Hence the probability of this
> sequence is (1/2)^12 = 0.00024.  So the p-value of the (one-sided)  
> test
> is 0.00024.  Hence the result is ``significant'' at the usual levels,
> and my vet friend is happy.
> I would very much appreciate comments on my reasoning.  Have I made  
> any
> goof-ups, missed any obvious pit-falls?  Gone down a wrong garden  
> path?
> Is there a better approach?
> Most importantly (!!!): Is there any literature in which this approach
> is spelled out?  (The journal in which she wishes to publish will  
> almost
> surely demand a citation.  They *won't* want to see the reasoning
> spelled out in the paper.)
> I would conjecture that this sort of scenario must arise reasonably
> often in medical statistics and the suggested approach (if it is  
> indeed
> valid and sensible) would be ``standard''.  It might even have a name!
> But I have no idea where to start looking, so I thought I'd ask this
> wonderfully learned list.
> Thanks for any input.
> 	cheers,
> 		Rolf Turner
> ######################################################################
> Attention:\ This e-mail message is privileged and confid... 
> {{dropped:9}}
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

Attention:\ This e-mail message is privileged and confid...{{dropped:9}}

More information about the R-help mailing list