dunnett.PB {pbANOVA}

PB multiple comparisons of the levels of factor A aginst the control group

Description

Using Parametric Bootstrap to simulate a distribution for the multiple comparisons of treatment groups against a control

Usage

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

Value

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

#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