1. stem bmi
Stem-and-leaf plot for bmi
bmi rounded to nearest multiple of .1
plot in units of .1
18* | 6
19* | 3
20* |
21* | 229
22* | 223
23* | 4445679
24* | 0455689
25* | 247799
26* | 1334677889
27* | 4457788
28* | 166
29* | 3566788
30* | 3379
31* | 2
32* | 12236
33* | 29
34* | 15
35* | 07
. lv bmi
# 71 bmi
------
M 36 | 26.8 | spread pseudosigma
F 18.5 | 24.5 27.125 29.75 | 5.25 3.978136
E 9.5 | 23.4 27.8 32.2 | 8.800001 3.881823
D 5 | 21.9 27.9 33.9 | 12 3.971424
C 3 | 21.2 27.85 34.5 | 13.3 3.731982
B 2 | 19.3 27.15 35 | 15.7 3.947206
A 1.5 | 18.95 27.15 35.35 | 16.4 3.839632
1 | 18.6 27.15 35.7 | 17.1 3.65948
| |
| | # below # above
inner fence | 16.625 37.625 | 0 0
outer fence | 8.75 45.5 | 0 0
. summ bmi,detail
bmi
------
Percentiles Smallest
1% 18.6 18.6
5% 21.2 19.3
10% 22.3 21.2 Obs 71
25% 24.5 21.2 Sum of Wgt. 71
50% 26.8 Mean 27.21127
Largest Std. Dev. 3.828896
75% 29.8 34.1
90% 32.3 34.5 Variance 14.66044
95% 34.1 35 Skewness .1755801
99% 35.7 35.7 Kurtosis 2.557761
.
No extreme points as there are no values beyond the fences. See the letter value command (lv bmi)
Most would agree that these graphs and summaries offer no indications of asymmetry. As bmi is a ratio, bmi often does display skewness, but this example is a sample from a very obese population and is certainly not like a ‘general’ population of bmi values.
Have a look at these ‘tests for normality’
. sktest bmi
Skewness/Kurtosis tests for Normality
------joint ------
Variable | Pr(Skewness) Pr(Kurtosis) adj chi2(2) Prob>chi2
------+------
bmi | 0.514 0.511 0.88 0.6441
. swilk bmi
Shapiro-Wilk W test for normal data
Variable | Obs W V z Prob>z
------+------
bmi | 71 0.98772 0.764 -0.585 0.72071
In Stata, look under distributional graphs, below is an example of a symmetry plot (symplot bmi)
Next is an example of a normal probability plot. (pnorm bmi)
2.
. list
+------+
| resp group |
|------|
1. | 0 0 |
2. | 0 0 |
3. | 0 0 |
4. | 0 0 |
5. | 1 1 |
|------|
6. | 1 1 |
7. | 1 1 |
8. | 3 1 |
+------+
. permute group `"regress resp group"' _b, reps(70)
command: regress resp group
statistics: b_group = _b[group]
b_cons = _b[_cons]
permute var: group
Monte Carlo permutation statistics Number of obs = 8
Replications = 70
------
T | T(obs) c n p=c/n SE(p) [95% Conf. Interval]
------+------
b_group | 1.5 2 70 0.0286 0.0199 .0034791 .0994289
b_cons | -2.22e-16 70 70 1.0000 0.0000 .9486662 1
------
Note: confidence intervals are with respect to p=c/n
Note: c = #{|T| >= |T(obs)|}
The example above illustrates an estimation of the p-value. By fluke, I got the ‘exact’ p-value.
a)
0000 1113 1 –1.5
0001 0113 12 –1.0 (12 is 4x3: 4C3 x 3C1) choose 3 of the 4 0s and 1 of the 3 1s
0003 0111 4 0.0
0011 0013 18 –0.5 (18 is 6x3: 4c2 x 3C2)
0013 0011 18 0.5
0111 0003 4 0.0
0113 0001 12 1.0
1113 0000 1 1.5
b) symmetric about 0 (note you should add 4+4=8 to determine the count for 0 i.e. the number of ways the diff can be zero comes from 2 different rearrangements)
diff / -1.5 / -1 / -0.5 / 0 / 0.5 / 1 / 1.5prob / 1/70 / 12/70 / 18/70 / 8/70 / 18/70 / 12/70 / 1/70
A worthwhile check is that 1+12+18+8+12+1 = 70
The probabilities in the probability distribution must sum to one.
c) 2-sided p-value is 2/70
d) Null hypothesis: 1 = 2 where 1 is the mean for group 1 and 2 is the mean for group 2. In ‘real’ studies, you would indicate mean ‘what’. I.e. mean weight, mean systolic blood pressure, mean HbA1c etc…
(It is incorrect to pose a null hypothesis such as X_1 = X_2 : we know from our data that X_1is not X_2 . Hypotheses must be specified in terms of population characteristics: means, medians, probabilities, RR, OR, RD and so on… )
At the 5% level of significance, there is evidence against the null hypothesis that 1 = 2 (p-value is 2/70 = .0286)
The number of decimal places for the p-value is usually dictated by the journal. Noting p < 0.05 is very old fashioned and actual values should be given whenever possible. The only exception is when the p-value is very small. In this case it is best to write (for example) p < 0.0001. Again, the number of decimal places for such inequalities is usually standardized by the journal.
3. anova red dose
Number of obs = 60 R-squared = 0.0948
Root MSE = 13.8329 Adj R-squared = 0.0631
Source | Partial SS df MS F Prob > F
------+------
Model | 1142.8 2 571.4 2.99 0.0584
|
dose | 1142.8 2 571.4 2.99 0.0584
|
Residual | 10906.85 57 191.348246
------+------
Total | 12049.65 59 204.231356
. matrix c1 = (0,1,-2,1)
. matrix c2 = (0,-1,0,1)
. matrix c=c1\c2
. test,test(c) mtest(noadj)
( 1) dose[1] - 2 dose[2] + dose[3] = 0
( 2) - dose[1] + dose[3] = 0
------
| F(df,57) df p
------+------
(1) | 5.09 1 0.0279 # <- notice t*t = F
(2) | 0.88 1 0.3526 # <- not used here
------+------
all | 2.99 2 0.0584
------
# unadjusted p-values
. disp 0.25-2*10.85+4.35
-17.1
. disp 13.8329*sqrt(6/20)
7.5765914
. disp (0.25-2*10.85+4.35)/(13.8329*sqrt(6/20))
-2.2569516
. disp -2.2569516*-2.2569516
5.0938305
. disp 2*ttail(57,2.2569516)
.02786124
.
. table dose lin
------
| lin
dose | -1 0 1
------+------
.03 | 20
.04 | 20
.05 | 20
------
. table dose qua
------
| qua
dose | -2 1
------+------
.03 | 20
.04 | 20
.05 | 20
------
. gen qc=qua/6
. gen lc=lin/2
. regr red lc qc
Source | SS df MS Number of obs = 60
------+------F( 2, 57) = 2.99
Model | 1142.8 2 571.4 Prob > F = 0.0584
Residual | 10906.85 57 191.348246 R-squared = 0.0948
------+------Adj R-squared = 0.0631
Total | 12049.65 59 204.231356 Root MSE = 13.833
------
red | Coef. Std. Err. t P>|t| [95% Conf. Interval]
------+------
lc | 4.1 4.374337 0.94 0.353 -4.659459 12.85946
qc | -17.1 7.576574 -2.26 0.028 -32.27183 -1.928172
_cons | 5.15 1.785816 2.88 0.006 1.573966 8.726034
------
1 is the mean reduction in blood pressure for those patients receiving a 3% dose
2 is the mean reduction for those receiving 4%
and 3 is the mean reduction for those receiving 5%.
These 3 means are points on the curve of mean reduction versus dose. This curve is the (unknown, population) mean dose response curve.
The null hypothesis is: the mean dose response curve is linear. OR the mean dose response curve is a straight line. OR the 3 means are collinear.
OR = 1 - 2*2 +3 = 0
The estimate of is g. (g is the latin equivalent of (gamma))
The p-value is 2*Prob(t572.256) = 0.0279
At the 5% level of significance, there is evidence that the mean dose response curve is nonlinear.
The next step in this type of analysis would probably be to provide a confidence interval for the maximum on the dose response curve assuming a parabola describes the form of this curve. We will address this later in the course.
It is generally regarded as best to stop any further statistical testing so as to not diminish the impact the nonlinearity detection.
Other confidence intervals could be offered though. Many were suggested by the class.
4. csi 7 4 2 8,exact
| Exposed Unexposed | Total
------+------+------
Cases | 7 4 | 11
Noncases | 2 8 | 10
------+------+------
Total | 9 12 | 21
| |
Risk | .7777778 .3333333 | .5238095
| |
| Point estimate | [95% Conf. Interval]
|------+------
Risk difference | .4444444 | .0637727 .8251162
Risk ratio | 2.333333 | .9745889 5.586401
Attr. frac. ex. | .5714286 | -.0260736 .8209939
Attr. frac. pop | .3636364 |
+------
1-sided Fisher's exact P = 0.0563
2-sided Fisher's exact P = 0.0805
. cci 7 4 2 8,exact
Proportion
| Exposed Unexposed | Total Exposed
------+------+------
Cases | 7 4 | 11 0.6364
Controls | 2 8 | 10 0.2000
------+------+------
Total | 9 12 | 21 0.4286
| |
| Point estimate | [95% Conf. Interval]
|------+------
Odds ratio | 7 | .7307428 90.81204 (exact)
Attr. frac. ex. | .8571429 | -.3684705 .9889882 (exact)
Attr. frac. pop | .5454545 |
+------
1-sided Fisher's exact P = 0.0563
2-sided Fisher's exact P = 0.0805
P(D|E) = Probability of No remission given mouse received 6-MP
P(D| not E) = Probability of No Remission given mouse received mGAG
Risk Difference = RD = P(D|E) – P(D| not E)
Risk Ratio = P(D|E)/P(D| not E)
Odds of No Remission given mouse received 6-MP = Odds(D|E) =
P(D|E)/P(not D |E)
Odds of No Remission given mouse received mGAG = Odds(D| not E)
P(D| not E)/P(not D |not E)
Odds Ratio = OR = Odds(D|E)/Odds(D| not E)
Estimate of RD = R^D = 7/9 – 4/12 = 4/9
95% CI for RD (0.0637, 0.8252) from Stata
Estimate of RR = R^R = (7/9) /(4/12) = 7/3
95% CI for RR (0.9745, 5.5865)
Estimate of OR = O^R = (7/2) /(4/8) = 7
95% CI for OR (0.7307, 90.8121)
Null hypothesis: RD = 0 OR RR=1 OR OR = 1 (terms are defined above)
The2test here is an approximate test. Since the exact test is available, it would be inappropriate to interpret the result of this test. Exact tests are always preferred over approximate tests as we are then able to interpret the probability without concern over cell sizes or obsolete ‘rules of thumb’.
Confidence intervals are helpful here as included above. It would be ill advised to try a stratified analysis here as the cells are already small. In principle though, assessment of potential confounders and/or potential modifiers would make sense for the next (much larger) study.
This little study assessed here is surely not much more than a pilot study.
5. See Sleuth page 107/108 question 28
. list
+------+
| pair ht1 ht2 diff |
|------|
1. | 1 188 139 49 |
2. | 2 96 163 -67 |
3. | 3 168 160 8 |
4. | 4 176 160 16 |
5. | 5 153 147 6 |
|------|
6. | 6 172 149 23 |
7. | 7 177 149 28 |
8. | 8 163 122 41 |
9. | 9 146 132 14 |
10. | 10 173 144 29 |
|------|
11. | 11 186 130 56 |
12. | 12 168 144 24 |
13. | 13 177 102 75 |
14. | 14 184 124 60 |
15. | 15 96 144 -48 |
+------+
. stem diff,lines(1)
Stem-and-leaf plot for diff
-6* | 7
-5* |
-4* | 8
-3* |
-2* |
-1* |
-0* |
0* | 68
1* | 46
2* | 3489
3* |
4* | 19
5* | 6
6* | 0
7* | 5
. lv diff
# 15 diff
------
M 8 | 24 | spread pseudosigma
F 4.5 | 11 28 45 | 34 27.9807
E 2.5 | -21 18.5 58 | 79 36.76205
D 1.5 | -57.5 5 67.5 | 125 43.64846
1 | -67 4 75 | 142 41.91955
| |
| | # below # above
inner fence | -40 96 | 2 0
outer fence | -91 147 | 0 0
. summ diff,detail
diff
------
Percentiles Smallest
1% -67 -67
5% -67 -48
10% -48 6 Obs 15
25% 8 8 Sum of Wgt. 15
50% 24 Mean 20.93333
Largest Std. Dev. 37.74438
75% 49 49
90% 60 56 Variance 1424.638
95% 75 60 Skewness -.9920929
99% 75 75 Kurtosis
a)See the actual description in the text. Many acceptable descriptions offered.
b)If d is the mean difference in height between cross fertilized and self fertilized plants then the null hypothesis is
H0: d = 0
c)The p-value = 0.0497. The 95% CI for d is (0.0311, 41.8355)
d)The distribution of the differences is normal. The differences are independent.
e)and f) Pairs 2 and 15 have differences (-48, –67) that are beyond the lower inner fence of –40.
g)If these 2 differences are removed, the average difference goes up, the sample standard deviation goes down, the t-statistic goes up as it is a ratio of something bigger divided by something smaller, the p-value goes down, the confidence intervals become narrower and the lower limit would be further from the null value (0)
. gen diff2 = diff
. replace diff2=. if diff2<-40
(2 real changes made, 2 to missing)
. ttest diff2=0
One-sample t test
------
Variable | Obs Mean Std. Err. Std. Dev. [95% Conf. Interval]
------+------
diff2 | 13 33 5.974304 21.54066 19.98311 46.01689
------
Degrees of freedom: 12
Ho: mean(diff2) = 0
Ha: mean < 0 Ha: mean != 0 Ha: mean > 0
t = 5.5237 t = 5.5237 t = 5.5237
P < t = 0.9999 P > |t| = 0.0001 P > t = 0.0001
. signtest diff=0
Sign test
sign | observed expected
------+------
positive | 13 7.5
negative | 2 7.5
zero | 0 0
------+------
all | 15 15
One-sided tests:
Ho: median of diff = 0 vs.
Ha: median of diff > 0
Pr(#positive >= 13) =
Binomial(n = 15, x >= 13, p = 0.5) = 0.0037
Ho: median of diff = 0 vs.
Ha: median of diff < 0
Pr(#negative >= 2) =
Binomial(n = 15, x >= 2, p = 0.5) = 0.9995
Two-sided test:
Ho: median of diff = 0 vs.
Ha: median of diff != 0
Pr(#positive >= 13 or #negative >= 13) =
min(1, 2*Binomial(n = 15, x >= 13, p = 0.5)) = 0.0074
The sign test has a p-value that is similar to that found with the t-test after removal of the 2 extreme points and is not disturbed by the 2 extreme points. Confidence intervals for the median can be obtained by other methods, but not directly from the sign test, per se.
. signrank diff=0
Wilcoxon signed-rank test
sign | obs sum ranks expected
------+------
positive | 13 96 60
negative | 2 24 60
zero | 0 0 0
------+------
all | 15 120 120
unadjusted variance 310.00
adjustment for ties 0.00
adjustment for zeros 0.00
------
adjusted variance 310.00
Ho: diff = 0
z = 2.045
Prob > |z| = 0.0409
The signrank test has a p-value that is very similar to the ‘original’ t-test.
h)The average difference would multiplied by 1/8. The sample standard deviation of differences would be multiplied by 1/8. All the letter values would be multiplied by 1/8. The t-statistic would then be unchanged and the p-value would be unchanged. The limits of the confidence interval for d would be multiplied by 1/8.
. gen diffinch=diff/8
. ttest diffinch=0
One-sample t test
------
Variable | Obs Mean Std. Err. Std. Dev. [95% Conf. Interval]
------+------
diffinch | 15 2.616667 1.218195 4.718047 .0038992 5.229434
------
Degrees of freedom: 14
Ho: mean(diffinch) = 0
Ha: mean < 0 Ha: mean != 0 Ha: mean > 0
t = 2.1480 t = 2.1480 t = 2.1480
P < t = 0.9751 P > |t| = 0.0497 P > t = 0.0249
. lv diffinch
# 15 diffinch
------
M 8 | 3 | spread pseudosigma
F 4.5 | 1.375 3.5 5.625 | 4.25 3.497588
E 2.5 | -2.625 2.3125 7.25 | 9.875 4.595256
D 1.5 | -7.1875 .625 8.4375 | 15.625 5.456057
1 | -8.375 .5 9.375 | 17.75 5.239943
| |
| | # below # above
inner fence | -5 12 | 2 0
outer fence | -11.375 18.375 | 0 0
. summ diffinch,detail
diffinch
------
Percentiles Smallest
1% -8.375 -8.375
5% -8.375 -6
10% -6 .75 Obs 15
25% 1 1 Sum of Wgt. 15
50% 3 Mean 2.616667
Largest Std. Dev. 4.718047
75% 6.125 6.125
90% 7.5 7 Variance 22.25997
95% 9.375 7.5 Skewness -.9920929
99% 9.375 9.375 Kurtosis 3.60574
6. table before
------
before | Freq.
------+------
1 | 1
2 | 5
3 | 15
4 | 20
5 | 29
6 | 38
7 | 30
8 | 28
9 | 13
10 | 8
11 | 3
12 | 2
13 | 1
14 | 1
25 | 1
30 | 1
33 | 1
35 | 1
120 | 1
200 | 1
------
. table after
------
after | Freq.
------+------
1 | 7
2 | 28
3 | 43
4 | 47
5 | 28
6 | 16
7 | 16
8 | 4
10 | 5
15 | 2
18 | 1
20 | 1
70 | 1
100 | 1
------
. lv before after
# 200 before
------
M 100.5 | 6 | spread pseudosigma
F 50.5 | 5 6.5 8 | 3 2.228226
E 25.5 | 4 6.5 9 | 5 2.178984
D 13 | 3 7 11 | 8 2.61741
C 7 | 3 8.5 14 | 11 2.997827
B 4 | 2 17.5 33 | 31 7.415788
A 2.5 | 2 39.75 77.5 | 75.5 16.43597
Z 1.5 | 1.5 80.75 160 | 158.5 31.41531
1 | 1 100.5 200 | 199 36.85351
| |
| | # below # above
inner fence | .5 12.5 | 0 8
outer fence | -4 17 | 0 6
# 200 after
------
M 100.5 | 4 | spread pseudosigma
F 50.5 | 3 4 5 | 2 1.485484
E 25.5 | 2 4.5 7 | 5 2.178984
D 13 | 2 5 8 | 6 1.963058
C 7 | 1 5.5 10 | 9 2.452767
B 4 | 1 9.5 18 | 17 4.066723
A 2.5 | 1 23 45 | 44 9.578576
Z 1.5 | 1 43 85 | 84 16.64912
1 | 1 50.5 100 | 99 18.33416
| |
| | # below # above
inner fence | 0 8 | 0 15
outer fence | -3 11 | 0 6
. lv lbefore lafter
# 200 lbefore
------
M 100.5 | 1.791759 | spread pseudosigma
F 50.5 | 1.609438 1.84444 2.079442 | .4700036 .3490914
E 25.5 | 1.386294 1.791759 2.197225 | .8109303 .3534009
D 13 | 1.098612 1.748254 2.397895 | 1.299283 .4250946
C 7 | 1.098612 1.868835 2.639057 | 1.540445 .4198171
B 4 | .6931472 2.094827 3.496508 | 2.80336 .670617
A 2.5 | .6931472 2.432284 4.17142 | 3.478273 .7572023
Z 1.5 | .3465736 2.694739 5.042905 | 4.696331 .9308309
1 | 0 2.649159 5.298317 | 5.298317 .9812139
| |
| | # below # above
inner fence | .9044325 2.784447 | 6 6
outer fence | .1994271 3.489452 | 1 4
# 200 lafter
------
M 100.5 | 1.386294 | spread pseudosigma
F 50.5 | 1.098612 1.354025 1.609438 | .5108256 .3794116
E 25.5 | .6931472 1.319529 1.94591 | 1.252763 .5459501
D 13 | .6931472 1.386294 2.079442 | 1.386294 .4535626
C 7 | 0 1.151293 2.302585 | 2.302585 .6275229
B 4 | 0 1.445186 2.890372 | 2.890372 .6914318
A 2.5 | 0 1.811057 3.622114 | 3.622114 .7885157
Z 1.5 | 0 2.213416 4.426833 | 4.426833 .8774153
1 | 0 2.302585 4.60517 | 4.60517 .8528476
| |
| | # below # above
inner fence | .3323739 2.375676 | 7 6
outer fence | -.4338646 3.141915 | 0 2
.
a)Huge outliers, clear skewness, difference in SDs, discrete data (days), average > median in both groups. A logarithmic transformation may reduce the skewness and bring the SDs closer together. Boxplots of log duration STILL show outliers. Must check these durations? Are they real? Analysis with and without them? Many options. None great. Data trouble, for sure.
b) The null hypothesis could be specified 2 ways:
In terms of the original units (days),
H0: the population median duration after the intervention =
the population median duration before the intervention
Or, in terms of the logarithmic units (log days)
H0: the mean log duration after = the mean log duration before
BUT not in terms of mean duration. It is clear that mean duration is not a meaningful measure of the centre of the duration distribution.
c) The p-value from Stata reads ‘= 0.0000’. This value should noted as p-value < 0.00005.
The 95% confidence interval for the difference in mean log duration is: (-0.3377, -0.5671). (this is before minus after)
The 95% confidence interval for the ratio of median durations is the exponent (with base e) which is: (0.5671, 0.7135) (this is before divided by after)
With 95% confidence, the median duration after intervention is at most 72% of the median duration before intervention. The median duration of stay in hospital divides the lengths of stay into the upper half and the lower half of the population distribution of durations indicating here a meaningful reduction in length of stay after intervention compared to before intervention.
However, data integrity issues have been identified in part a) and a thorough cleaning of this data seems to be in order before a ‘final report’ is distributed.