# PB multiple comparisons of the levels of factor A aginst the control group Using Parametric Bootstrap to simulate a distribution for the multiple comparisons of treatment groups against a control ```r dunnett.PB(L, ns, means, s2, alpha) ``` ## Arguments - `L`: Number of simulated values for the distribution - `ns`: sample size for each group - `means`: sample mean for each group - `s2`: sample variance for each group - `alpha`: significant level ## Returns D.crit: The (1 - alpha) percentile of the simulated distribution result: The differences, confidence intervals for the difference, and p-values for comparisons of each factor level vs. the control. ## Examples ```r #This one gets a different result between the PB method and the traditional Dunnett's test. #Constant variance assumption appears violated on residual plots. #The breusch pagan test shows close to violating (p=0.0596) while levene's test (using #median) #does not show a violation. Data is mildly unbalanced. #Traditional Dunnett's test says group 4-6 are different from "control", while the PB method # only identifies group 5 and 6. #Group 4 has a larger variance than the others. The pooled variance/MSE could be too small for #this group and lead to arificially large test statistic. #MSE is 0.0004274 and the sample variance of group 4 is 0.00169. #The authors of the paper do not claim significance, they report the means and state #that there is a delineation between 30 and 40 feet. This seems true when looking at the #means, but the measurements at 40 feet do have a larger variance than those at other depths. #Practical interpretation: #If your goal was to get the most iron rich water from as shallow depth as possible, #knowing that the surface (control?) was not rich enough, and you decided to go 40 feet # deep, you may still get water that didn't have enough iron content for your purpose. #Suppose you wanted at least 0.1 content; based on the means you might use 40 feet, but # from their data, the measurements were: #0.098, 0.074, 0.154 #So 2/3 of these samples wouldn't contain enough iron for your purpose. library(DescTools) ##Dunnett's test library(lmtest) #BP test for constant variance library(dplyr) #data manipulation library(MASS) fedata$depth <- factor(fedata$depth) femod <- lm(Y~depth, data=fedata) plot(femod$fit, rstandard(femod), main="Fitted-Residual Plot, One-Way ANOVA Model", sub="Iron Data") #appears to violate equal variance assumption bptest(femod) #close to violation #what about normality? qqnorm(rstandard(femod), main="Normal QQ-Plot, Standardized Residuals", sub="One-Way ANOVA Model, Iron Data") shapiro.test(rstandard(femod)) #does not violate fe.sums <-fedata %>% group_by(depth) %>% summarise(means=mean(Y), vars = var(Y), sd=sd(Y), ns=length(depth)) fe.sums #summarySE(fedata, "Y", "depth") from Rmisc pkg does same thing as fe.sums above pbd.fe <- dunnett.PB(L=5000, ns=fe.sums$ns, means=fe.sums$means, s2=fe.sums$vars, alpha=0.05) pbd.fe$result pbd.fe$D.crit ##grp 5 and 6 sig diff from group 1 DunnettTest(Y~depth, data=fedata) ##Dunnett's test also says group 4 is different. ##group 4 has a larger variance than the others. ##pooled variance/MSE could be too small for this group ##and lead to arificially large test statistic. anova(femod) #MSE is 0.0004274 #sample variance of group 4 is 0.00169 #Use traditional Dunnett's test on transformed data for comparisons #attempt log transformation fedata$logY <-log(fedata$Y) felogmod <- lm(logY~depth, data=fedata) bptest(felogmod) #still violates #LeveneTest(logY~depth, data=fedata) plot(felogmod$fit, rstandard(felogmod)) shapiro.test(rstandard(felogmod)) #still for normality #square root transformation? fedata$srY <- sqrt(fedata$Y) fesrmod <- lm(srY~depth, data=fedata) bptest(fesrmod) #still violates, worse plot(fesrmod$fit, rstandard(fesrmod)) boxcox(femod) #lambda=-0.2 fedata$bcY <- with(fedata, (Y^(-0.2) - 1)/-0.2) febcmod <- lm(bcY~depth, data=fedata) #gives same F-statistic as lm(Y^(-0.2) ~ depth, data=fedata), bptest(febcmod) #still violates plot(febcmod$fit, rstandard(febcmod)) ##mg/L is a proportion, try arcsin fedata$asY <- asin(sqrt(fedata$Y)) feasmod <-lm(asY~depth, data=fedata) bptest(feasmod) #still close to violating plot(felogmod$fit, rstandard(felogmod), main="Log Transform.", xlab="Fitted Values", ylab="Standardized Residuals") plot(febcmod$fit, rstandard(febcmod), main="Box-Cox Transform.", xlab="Fitted Values", ylab="Standardized Residuals") plot(feasmod$fit, rstandard(feasmod), main="Arcsin(Sq. Root) Transform.", xlab="Fitted Values", ylab="Standardized Residuals") #P-values of BP test are similar for log and box-cox, plots look a little better ##log transform may be considered simpler, so try that anova(felogmod) DunnettTest(logY~depth, data=fedata) #still identifies 40 feet and above DunnettTest(bcY~depth, data=fedata) #still identifies 40 feet and above shapiro.test(rstandard(febcmod)) # normality still ok #W = 0.95535, p-value = 0.4556 ```

dunnett.PB function