4.3.1 Estimating the Width of a Room Revised
The unconditional analysis of the room width estimated by two groups ofstudents in Chapter˜3 led to the conclusion that the estimates in metres areslightly larger than the estimates in feet. Here, we reanalyse these data in aconditional framework. First, we convert metres into feet and store the vectorof observations in a variable y:
R> data("roomwidth", package = "HSAUR2")R> convert <- ifelse(roomwidth$unit == "feet", 1, 3.28)R> feet <- roomwidth$unit == "feet"R> metre <- !feetR> y <- roomwidth$width * convert
The test statistic is simply the difference in means
R> T <- mean(y[feet]) - mean(y[metre])R> T
In order to approximate the conditional distribution of the test statistic Twe compute 9999 test statistics for shuffled y values. A permutation of the yvector can be obtained from the sample function.
R> meandiffs <- double(9999)R> for (i in 1:length(meandiffs)) {+
meandiffs[i] <- mean(sy[feet]) - mean(sy[metre])
The distribution of the test statistic T under the null hypothesis of indepen-dence of room width estimates and groups is depicted in FigureNow, thevalue of the test statistic T for the original unshuffled data can be compared
R> hist(meandiffs)R> abline(v = T, lty = 2)R> abline(v = -T, lty = 2)
Histogram of meandiffs
An approximation for the conditional distribution of the difference ofmean roomwidth estimates in the feet and metres group under the nullhypothesis. The vertical lines show the negative and positive absolutevalue of the test statistic T obtained from the original data.
with the distribution of T under the null hypothesis (the vertical lines in Fig-ure. The p-value, i.e., the proportion of test statistics T larger than 8.859or smaller than -8.859, is
R> greater <- abs(meandiffs) > abs(T)R> mean(greater)
R> binom.test(sum(greater), length(greater))$conf.int
Note that the approximated conditional p-value is roughly the same as thep-value reported by the t-test in Chapter˜3.
R> library("coin")R> independence_test(y ~ unit, data = roomwidth,+
R output of the exact permutation test applied to the roomwidth data.
R> wilcox_test(y ~ unit, data = roomwidth,+
Exact Wilcoxon Mann-Whitney Rank Sum Testalternative hypothesis: true mu is not equal to 0
R output of the exact conditional Wilcoxon rank sum test applied tothe roomwidth data.
Here we are interested in the comparison of two groups of patients, where onegroup received a placebo and the other one Misoprostol. In the trials shownhere, the response variable is measured on an ordered scale – see Table˜??. Data from four clinical studies are available and thus the observations arenaturally grouped together. From the data.frame Lanza we can construct athree-way table as follows:
R> data("Lanza", package = "HSAUR2")R> xtabs(~ treatment + classification + study, data = Lanza)
R> data("suicides", package = "HSAUR2")R> fisher.test(suicides)
Fisher's Exact Test for Count Dataalternative hypothesis: true odds ratio is not equal to 1
R output of Fisher’s exact test for the suicides data.
For the first study, the null hypothesis of independence of treatment and
gastrointestinal damage, i.e., of no treatment effect of Misoprostol, is testedby
R> library("coin")R> cmh_test(classification ~ treatment, data = Lanza,+
scores = list(classification = c(0, 1, 6, 17, 30)),
Asymptotic Linear-by-Linear Association Testclassification (ordered) by treatment (Misoprostol, Placebo)chi-squared = 28.8, df = 1, p-value = 7.83e-08
and, by default, the conditional distribution is approximated by the corre-sponding limiting distribution. The p-value indicates a strong treatment effect. For the second study, the asymptotic p-value is a little bit larger:
R> cmh_test(classification ~ treatment, data = Lanza,+
scores = list(classification = c(0, 1, 6, 17, 30)),
Asymptotic Linear-by-Linear Association Testclassification (ordered) by treatment (Misoprostol, Placebo)chi-squared = 12.1, df = 1, p-value = 0.000514
and we make sure that the implied decision is correct by calculating a confi-dence interval for the exact p-value:
R> p <- cmh_test(classification ~ treatment, data = Lanza,+
scores = list(classification = c(0, 1, 6, 17, 30)),
subset = Lanza$study == "II", distribution =
The third and fourth study indicate a strong treatment effect as well:
R> cmh_test(classification ~ treatment, data = Lanza,+
scores = list(classification = c(0, 1, 6, 17, 30)),
Asymptotic Linear-by-Linear Association Testclassification (ordered) by treatment (Misoprostol, Placebo)chi-squared = 28.2, df = 1, p-value = 1.118e-07
R> cmh_test(classification ~ treatment, data = Lanza,+
scores = list(classification = c(0, 1, 6, 17, 30)),
Asymptotic Linear-by-Linear Association Testclassification (ordered) by treatment (Misoprostol, Placebo)chi-squared = 15.7, df = 1, p-value = 7.262e-05
At the end, a separate analysis for each study is unsatisfactory. Because the
design of the four studies is the same, we can use study as a block variableand perform a global linear-association test investigating the treatment effectof Misoprostol in all four studies. The block variable can be incorporated intothe formula by the | symbol.
R> cmh_test(classification ~ treatment | study, data = Lanza,+
scores = list(classification = c(0, 1, 6, 17, 30)))
Asymptotic Linear-by-Linear Association Testchi-squared = 83.6, df = 1, p-value < 2.2e-16
Based on this result, a strong treatment effect can be established.
In this example, the medical doctor (MD) and the research assistant (RA)assessed the number of anomalies (0, 1, 2 or 3) for each of 395 babies:
R> anomalies <- c(235, 23, 3, 0, 41, 35, 8, 0,+
R> anomalies <- as.table(matrix(anomalies,+
ncol = 4, dimnames = list(MD = 0:3, RA = 0:3)))
We are interested in testing whether the number of anomalies assessed by themedical doctor differs structurally from the number reported by the researchassistant. Because we compare paired observations, i.e., one pair of measure-ments for each newborn, a test of marginal homogeneity (a generalisation ofMcNemar’s test, Chapter˜3) needs to be applied:
Asymptotic Marginal-Homogeneity Testchi-squared = 21.2, df = 3, p-value = 9.446e-05
The p-value indicates a deviation from the null hypothesis. However, the levelsof the response are not treated as ordered. Similar to the analysis of thegastrointestinal damage data above, we can take this information into accountby the definition of an appropriate score. Here, the number of anomalies is anatural choice:
R> mh_test(anomalies, scores = list(response = c(0, 1, 2, 3)))
Asymptotic Marginal-Homogeneity Test for Ordered Datachi-squared = 21, df = 1, p-value = 4.545e-06
In our case, one can conclude that the assessment of the number of anomaliesdiffers between the medical doctor and the research assistant.
Uncle Louis’ Guide to Case Write-Ups Most of this advice, though simple and logical is based largely on the comments I have found myself giving to students over the past few years about their assignment write-ups. Before one starts the assignment, one ought to consider some basics rules of engagement. Being that we are a Commonwealth Country, and Queen Elizabeth II of England, is also