lab6_bagga
Amneet Bagga
2023-03-21
Task 1
# Download the data
d1=read.csv("http://faraway.neu.edu/biostats/lab6_dataset1.csv")
d1## drug1 drug2 score
## 1 d1.ctrl d2.ctrl 130.783880
## 2 d1.ctrl d2.ctrl 294.018686
## 3 d1.ctrl d2.ctrl 106.098869
## 4 d1.ctrl d2.ctrl 1206.261037
## 5 d1.ctrl d2.ctrl 340.191186
## 6 d1.ctrl d2.ctrl 107.719607
## 7 d1.ctrl d2.ctrl 398.389055
## 8 d1.ctrl d2.ctrl 512.000041
## 9 d1.ctrl d2.ctrl 435.189406
## 10 d1.ctrl d2.ctrl 180.298104
## 11 d1.ctrl d2.ctrl 1109.629182
## 12 d1.ctrl d2.ctrl 361.348634
## 13 d1.ctrl d2.ctrl 131.467467
## 14 d1.ctrl d2.ctrl 26.717001
## 15 d1.ctrl d2.ctrl 753.652147
## 16 d1.ctrl d2.ctrl 233.940402
## 17 d1.ctrl d2.ctrl 240.762203
## 18 d1.ctrl d2.ctrl 628.814443
## 19 d1.ctrl d2.ctrl 556.251870
## 20 d1.ctrl d2.ctrl 443.146902
## 21 d1.trt d2.ctrl 1508.660382
## 22 d1.trt d2.ctrl 1315.716023
## 23 d1.trt d2.ctrl 648.437085
## 24 d1.trt d2.ctrl 82.322816
## 25 d1.trt d2.ctrl 1118.591684
## 26 d1.trt d2.ctrl 568.994780
## 27 d1.trt d2.ctrl 515.019361
## 28 d1.trt d2.ctrl 138.275437
## 29 d1.trt d2.ctrl 373.101293
## 30 d1.trt d2.ctrl 914.101449
## 31 d1.trt d2.ctrl 2341.810323
## 32 d1.trt d2.ctrl 543.055906
## 33 d1.trt d2.ctrl 886.846233
## 34 d1.trt d2.ctrl 570.318490
## 35 d1.trt d2.ctrl 151.857175
## 36 d1.trt d2.ctrl 397.424682
## 37 d1.trt d2.ctrl 405.738980
## 38 d1.trt d2.ctrl 567.185609
## 39 d1.trt d2.ctrl 1808.088289
## 40 d1.trt d2.ctrl 1291.004335
## 41 d1.ctrl d2.trt 110.088638
## 42 d1.ctrl d2.trt 104.374249
## 43 d1.ctrl d2.trt 184.597545
## 44 d1.ctrl d2.trt 169.694187
## 45 d1.ctrl d2.trt 80.378488
## 46 d1.ctrl d2.trt 79.479800
## 47 d1.ctrl d2.trt 151.222033
## 48 d1.ctrl d2.trt 192.697144
## 49 d1.ctrl d2.trt 113.589636
## 50 d1.ctrl d2.trt 206.162488
## 51 d1.ctrl d2.trt 154.294564
## 52 d1.ctrl d2.trt 84.165421
## 53 d1.ctrl d2.trt 149.108140
## 54 d1.ctrl d2.trt 61.706060
## 55 d1.ctrl d2.trt 287.095301
## 56 d1.ctrl d2.trt 398.712234
## 57 d1.ctrl d2.trt 97.481938
## 58 d1.ctrl d2.trt 64.943602
## 59 d1.ctrl d2.trt 171.028768
## 60 d1.ctrl d2.trt 112.052469
## 61 d1.trt d2.trt 365.628489
## 62 d1.trt d2.trt 31.841167
## 63 d1.trt d2.trt 66.005585
## 64 d1.trt d2.trt 34.055861
## 65 d1.trt d2.trt 15.748211
## 66 d1.trt d2.trt 39.996514
## 67 d1.trt d2.trt 5.446871
## 68 d1.trt d2.trt 143.388089
## 69 d1.trt d2.trt 38.600041
## 70 d1.trt d2.trt 290.792999
## 71 d1.trt d2.trt 53.277256
## 72 d1.trt d2.trt 16.281892
## 73 d1.trt d2.trt 60.991003
## 74 d1.trt d2.trt 13.012395
## 75 d1.trt d2.trt 9.453326
## 76 d1.trt d2.trt 44.320452
## 77 d1.trt d2.trt 21.257465
## 78 d1.trt d2.trt 33.152076
## 79 d1.trt d2.trt 35.671117
## 80 d1.trt d2.trt 18.3655951.1
#two-way ANOVA
d1_anova2w<- aov(score ~ drug1 * drug2, d1)
summary(d1_anova2w)## Df Sum Sq Mean Sq F value Pr(>F)
## drug1 1 498377 498377 4.224 0.04329 *
## drug2 1 5016549 5016549 42.518 6.92e-09 ***
## drug1:drug2 1 1148512 1148512 9.734 0.00256 **
## Residuals 76 8966901 117986
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1d1_anova2wlog10<- aov(log10(score) ~ drug1 * drug2, d1)
summary(d1_anova2wlog10)## Df Sum Sq Mean Sq F value Pr(>F)
## drug1 1 0.316 0.316 2.283 0.135
## drug2 1 12.216 12.216 88.344 2.30e-14 ***
## drug1:drug2 1 3.714 3.714 26.861 1.75e-06 ***
## Residuals 76 10.509 0.138
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Transforming the data does change the model fit, and significantly.
1.2
d1aovresid = resid(d1_anova2w)
d1aovlog10resid = resid(d1_anova2wlog10)
par(mfrow=c(1,2))
qqnorm(d1aovresid)
qqline(d1aovresid, col = "red")
qqnorm(d1aovlog10resid)
qqline(d1aovlog10resid, col = "red")
The log-transofrmed data meets the assumptions and should be used for
these data.
1.3
Drug 1 alone does not have a significant effect, but drug 2 does, and both drugs in combination do as well.
1.4
mod.tukey = TukeyHSD(d1_anova2wlog10)
plot(mod.tukey)
#### 1.5 There is a significant difference in the difference of the mean
levels in viral loads between almost every group, with the exception of
the control and drug 1 treatment group. It really does seem as if drug 2
preforms the best here, as it has a lower viral load level in comparison
to most of the other treatments.
1.6
require(multcompView, quietly=TRUE)
#generate the means and standard errors using tapply()
means <- tapply(d1$score, list(d1$drug1, d1$drug2), mean)
stderrs <- tapply(d1$score, list(d1$drug1, d1$drug2), FUN = function(x) sd(x)/sqrt(length(x)))
hiv_labels<- multcompLetters(mod.tukey$`drug1:drug2`[, 'p adj'])$Letters
hiv_labels<- hiv_labels[order(names(hiv_labels))]
print(hiv_labels)## d1.ctrl:d2.ctrl d1.ctrl:d2.trt d1.trt:d2.ctrl d1.trt:d2.trt
## "a" "b" "a" "c"#generate the barplot
hiv_bp<- barplot(means, beside = T, ylab = "Viral Load", xlab = "Drug 2 Treatment",
ylim = c(0, max(means) + 30), col = c("pink", "blue"))
arrows(x0 = hiv_bp, x1 = hiv_bp, y0 = means - stderrs, y1 = means + stderrs,
angle = 90, len = 0.1, code = 3)
text(x = hiv_bp, y = means + stderrs, pos = 3, labels = hiv_labels)
legend(x = "topleft", legend = c("drug1.ctrl", "drug1.trt"),
fill = c("pink", "blue"), title = "Drug 1 Treatment Group")
1.7
interaction.plot(d1$drug1, d1$drug2, d1$score)
I prefer the bar plot. It was easier for me to understand by looking really quickly.
1.8
A cocktail of both drugs is clearly the best option. Drug 1 seems to have no/and adverse effect, and drug 2 seems to be effective, but both drugs in combination seems to be the most effective. If we were to use one drug, however, I would very heavily lean drug 2.
Task 2
# Download the data
d2=read.csv("http://faraway.neu.edu/biostats/lab6_dataset2.csv")
d2## drug1 drug2 score sex
## 1 d1.ctrl d2.ctrl 99.791744 male
## 2 d1.ctrl d2.ctrl 294.373446 male
## 3 d1.ctrl d2.ctrl 1197.325180 male
## 4 d1.ctrl d2.ctrl 79.013943 male
## 5 d1.ctrl d2.ctrl 225.822263 male
## 6 d1.ctrl d2.ctrl 279.337376 male
## 7 d1.ctrl d2.ctrl 496.684358 male
## 8 d1.ctrl d2.ctrl 192.539625 male
## 9 d1.ctrl d2.ctrl 1780.187432 male
## 10 d1.ctrl d2.ctrl 212.983136 male
## 11 d1.ctrl d2.ctrl 371.537853 male
## 12 d1.ctrl d2.ctrl 653.114708 male
## 13 d1.ctrl d2.ctrl 165.224416 male
## 14 d1.ctrl d2.ctrl 86.516143 male
## 15 d1.ctrl d2.ctrl 1454.225827 male
## 16 d1.ctrl d2.ctrl 24.262475 male
## 17 d1.ctrl d2.ctrl 589.105085 male
## 18 d1.ctrl d2.ctrl 253.612299 male
## 19 d1.ctrl d2.ctrl 673.729498 male
## 20 d1.ctrl d2.ctrl 377.007528 male
## 21 d1.trt d2.ctrl 4869.853849 male
## 22 d1.trt d2.ctrl 181.285689 male
## 23 d1.trt d2.ctrl 2950.229374 male
## 24 d1.trt d2.ctrl 4249.903867 male
## 25 d1.trt d2.ctrl 604.824163 male
## 26 d1.trt d2.ctrl 51.846821 male
## 27 d1.trt d2.ctrl 969.942996 male
## 28 d1.trt d2.ctrl 331.438354 male
## 29 d1.trt d2.ctrl 1329.028190 male
## 30 d1.trt d2.ctrl 804.030103 male
## 31 d1.trt d2.ctrl 1260.090224 male
## 32 d1.trt d2.ctrl 827.956321 male
## 33 d1.trt d2.ctrl 1765.456107 male
## 34 d1.trt d2.ctrl 452.977420 male
## 35 d1.trt d2.ctrl 276.808166 male
## 36 d1.trt d2.ctrl 331.736010 male
## 37 d1.trt d2.ctrl 107.127554 male
## 38 d1.trt d2.ctrl 244.060347 male
## 39 d1.trt d2.ctrl 344.101986 male
## 40 d1.trt d2.ctrl 470.354861 male
## 41 d1.ctrl d2.trt 96.529462 male
## 42 d1.ctrl d2.trt 37.507394 male
## 43 d1.ctrl d2.trt 73.330525 male
## 44 d1.ctrl d2.trt 380.744475 male
## 45 d1.ctrl d2.trt 176.530975 male
## 46 d1.ctrl d2.trt 401.236994 male
## 47 d1.ctrl d2.trt 101.160642 male
## 48 d1.ctrl d2.trt 115.064571 male
## 49 d1.ctrl d2.trt 108.799107 male
## 50 d1.ctrl d2.trt 59.189215 male
## 51 d1.ctrl d2.trt 73.481062 male
## 52 d1.ctrl d2.trt 419.801010 male
## 53 d1.ctrl d2.trt 86.717158 male
## 54 d1.ctrl d2.trt 261.237331 male
## 55 d1.ctrl d2.trt 64.809774 male
## 56 d1.ctrl d2.trt 37.355235 male
## 57 d1.ctrl d2.trt 100.104771 male
## 58 d1.ctrl d2.trt 213.048005 male
## 59 d1.ctrl d2.trt 240.696759 male
## 60 d1.ctrl d2.trt 331.282425 male
## 61 d1.trt d2.trt 5.538689 male
## 62 d1.trt d2.trt 252.457399 male
## 63 d1.trt d2.trt 16.393021 male
## 64 d1.trt d2.trt 38.790088 male
## 65 d1.trt d2.trt 54.939622 male
## 66 d1.trt d2.trt 14.585165 male
## 67 d1.trt d2.trt 4.486859 male
## 68 d1.trt d2.trt 20.505793 male
## 69 d1.trt d2.trt 36.023803 male
## 70 d1.trt d2.trt 13.524642 male
## 71 d1.trt d2.trt 13.180314 male
## 72 d1.trt d2.trt 46.083248 male
## 73 d1.trt d2.trt 28.741417 male
## 74 d1.trt d2.trt 51.154362 male
## 75 d1.trt d2.trt 31.383346 male
## 76 d1.trt d2.trt 13.368345 male
## 77 d1.trt d2.trt 121.937941 male
## 78 d1.trt d2.trt 71.649758 male
## 79 d1.trt d2.trt 94.871714 male
## 80 d1.trt d2.trt 8.084605 male
## 81 d1.ctrl d2.ctrl 662.476172 female
## 82 d1.ctrl d2.ctrl 44.890900 female
## 83 d1.ctrl d2.ctrl 143.542026 female
## 84 d1.ctrl d2.ctrl 62.036973 female
## 85 d1.ctrl d2.ctrl 26.898761 female
## 86 d1.ctrl d2.ctrl 1513.412810 female
## 87 d1.ctrl d2.ctrl 127.307649 female
## 88 d1.ctrl d2.ctrl 184.070490 female
## 89 d1.ctrl d2.ctrl 166.176487 female
## 90 d1.ctrl d2.ctrl 360.212802 female
## 91 d1.ctrl d2.ctrl 1212.440866 female
## 92 d1.ctrl d2.ctrl 1314.425485 female
## 93 d1.ctrl d2.ctrl 74.917957 female
## 94 d1.ctrl d2.ctrl 62.899785 female
## 95 d1.ctrl d2.ctrl 53.910712 female
## 96 d1.ctrl d2.ctrl 69.888080 female
## 97 d1.ctrl d2.ctrl 1736.031564 female
## 98 d1.ctrl d2.ctrl 246.569986 female
## 99 d1.ctrl d2.ctrl 105.360184 female
## 100 d1.ctrl d2.ctrl 134.134080 female
## 101 d1.trt d2.ctrl 291.330804 female
## 102 d1.trt d2.ctrl 129.101360 female
## 103 d1.trt d2.ctrl 72.655422 female
## 104 d1.trt d2.ctrl 47.010449 female
## 105 d1.trt d2.ctrl 42.005547 female
## 106 d1.trt d2.ctrl 771.271385 female
## 107 d1.trt d2.ctrl 254.657646 female
## 108 d1.trt d2.ctrl 741.508864 female
## 109 d1.trt d2.ctrl 65.276130 female
## 110 d1.trt d2.ctrl 70.046939 female
## 111 d1.trt d2.ctrl 35.609747 female
## 112 d1.trt d2.ctrl 77.438252 female
## 113 d1.trt d2.ctrl 159.470582 female
## 114 d1.trt d2.ctrl 387.199406 female
## 115 d1.trt d2.ctrl 174.891993 female
## 116 d1.trt d2.ctrl 156.958288 female
## 117 d1.trt d2.ctrl 340.683421 female
## 118 d1.trt d2.ctrl 313.292371 female
## 119 d1.trt d2.ctrl 110.675007 female
## 120 d1.trt d2.ctrl 45.453205 female
## 121 d1.ctrl d2.trt 94.134056 female
## 122 d1.ctrl d2.trt 48.585010 female
## 123 d1.ctrl d2.trt 124.758386 female
## 124 d1.ctrl d2.trt 34.522643 female
## 125 d1.ctrl d2.trt 23.890055 female
## 126 d1.ctrl d2.trt 61.716278 female
## 127 d1.ctrl d2.trt 29.911661 female
## 128 d1.ctrl d2.trt 65.580826 female
## 129 d1.ctrl d2.trt 15.888667 female
## 130 d1.ctrl d2.trt 15.713842 female
## 131 d1.ctrl d2.trt 67.618560 female
## 132 d1.ctrl d2.trt 54.520562 female
## 133 d1.ctrl d2.trt 75.387436 female
## 134 d1.ctrl d2.trt 13.333140 female
## 135 d1.ctrl d2.trt 92.531346 female
## 136 d1.ctrl d2.trt 91.862861 female
## 137 d1.ctrl d2.trt 83.033239 female
## 138 d1.ctrl d2.trt 71.653798 female
## 139 d1.ctrl d2.trt 158.545895 female
## 140 d1.ctrl d2.trt 18.684870 female
## 141 d1.trt d2.trt 41.136514 female
## 142 d1.trt d2.trt 111.021512 female
## 143 d1.trt d2.trt 32.885057 female
## 144 d1.trt d2.trt 80.282402 female
## 145 d1.trt d2.trt 155.437414 female
## 146 d1.trt d2.trt 38.333605 female
## 147 d1.trt d2.trt 142.181724 female
## 148 d1.trt d2.trt 127.728796 female
## 149 d1.trt d2.trt 32.893354 female
## 150 d1.trt d2.trt 27.183775 female
## 151 d1.trt d2.trt 195.564576 female
## 152 d1.trt d2.trt 62.213799 female
## 153 d1.trt d2.trt 151.743528 female
## 154 d1.trt d2.trt 31.676700 female
## 155 d1.trt d2.trt 264.355208 female
## 156 d1.trt d2.trt 19.238038 female
## 157 d1.trt d2.trt 158.434407 female
## 158 d1.trt d2.trt 117.237591 female
## 159 d1.trt d2.trt 96.350567 female
## 160 d1.trt d2.trt 143.670529 female2.1
hiv.aov <- aov(d2$score ~ d2$drug1 * d2$drug2 * d2$sex)
summary(hiv.aov)## Df Sum Sq Mean Sq F value Pr(>F)
## d2$drug1 1 328367 328367 1.063 0.30422
## d2$drug2 1 8523822 8523822 27.588 4.99e-07 ***
## d2$sex 1 2597480 2597480 8.407 0.00429 **
## d2$drug1:d2$drug2 1 696277 696277 2.254 0.13539
## d2$drug1:d2$sex 1 1173594 1173594 3.798 0.05314 .
## d2$drug2:d2$sex 1 2092244 2092244 6.772 0.01018 *
## d2$drug1:d2$drug2:d2$sex 1 2540248 2540248 8.222 0.00473 **
## Residuals 152 46963816 308972
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 12.2
par(mfrow=c(1,2))
interaction.plot(d2$sex, d2$drug1, d2$score)
interaction.plot(d2$sex, d2$drug2, d2$score)
2.3
Drug 2 still seems the most effective, but soemthing really interesting happens when looking at drug 1. It performs effectively on females, but seems to have an adverse effect on men.
2.4
# Compute mean and standard error for each group
means <- aggregate(log10(d2$score) ~ d2$drug1 + d2$drug2 + d2$sex, FUN=mean, data=d2)
stderrs <- aggregate(log10(d2$score) ~ d2$drug1 + d2$drug2 + d2$sex, FUN = function(x) sd(x)/sqrt(length(x)), data=d2)
# Perform multiple comparisons
d2_anova2wlog10<- aov(log10(score) ~ drug1 * drug2 * sex, d2)
mod.aov.tukey=TukeyHSD(d2_anova2wlog10)
# Generate labels for each bar
labels=multcompLetters(mod.aov.tukey$`drug1:drug2:sex`[, "p adj"])$Letters
labels=labels[order(names(labels))]
means## d2$drug1 d2$drug2 d2$sex log10(d2$score)
## 1 d1.ctrl d2.ctrl female 2.270652
## 2 d1.trt d2.ctrl female 2.147756
## 3 d1.ctrl d2.trt female 1.694872
## 4 d1.trt d2.trt female 1.897523
## 5 d1.ctrl d2.ctrl male 2.473507
## 6 d1.trt d2.ctrl male 2.779362
## 7 d1.ctrl d2.trt male 2.109435
## 8 d1.trt d2.trt male 1.445205# Plot bars
bp=barplot(means$`log10(d2$score)`, beside = T, ylab="Viral Load", col=c("white", "lightgray", "darkgray", "black"))
arrows(x0 = bp, x1 = bp, y0 = means$`log10(d2$score)` - stderrs$`log10(d2$score)`, y1 = means$`log10(d2$score)` + stderrs$`log10(d2$score)`, angle = 90, len = 0.1, code = 3)
# Add labels for Males and Females
text(x=colMeans(bp)[1], y=-0.5, label='females males', xpd=NA)
box()
# Add bar labels based on multiple comparisons
text(x=bp, y=means[,4]+stderrs[,4], pos=3, lab=labels, xpd=NA)
# Create appropriate labels for the legend to describe the treatments
legend(x="topright", legend=c("d1.ctrl:d2.ctrl", "d1.trt:d2.ctrl", "d1.ctrl:d2.trt", "d1.trt:d2.trt"),
fill=c("white", "lightgray", "darkgray", "black"))
#### 2.5
The optimal treatment for male sand femlaes would be both drugs together. If we were to decide a different treamtment based on sex, the men would get a drug cocktail, but women would get just drug 2.