Normality, Variance and T-Tests
Homework assignment — Ruben Melikyan
This is a grades assignment
T1
Download and import the “Sprays.csv” file from Moodle. This dataset represents the counts of insects in agricultural experimental units treated with different insecticides.
All hypothesis tests below use a significance level of 0.05.
Sprays <- read.csv("Sprays.csv", row.names = 1)
Sprays$spray <- factor(Sprays$spray)
stopifnot(!anyNA(Sprays), is.numeric(Sprays$count))
head(Sprays)table(Sprays$spray)##
## A B C D E F
## 12 12 12 12 12 12There are 72 observations, with 12 for each spray.
T2
Make a Q-Q plot of the count variable from the above dataset. Make sure to add a 1:1 line (in red) (0.5 pt). Does it look like a normal distribution? Discuss (0.25 pt)
I standardize the counts so I can use a red 1:1 line. This does not change the shape of the data.
count_z <- as.numeric(scale(Sprays$count))
qqnorm(count_z, main = "Normal Q-Q plot of insect counts",
xlab = "Theoretical standard normal quantiles",
ylab = "Observed standardized insect counts", pch = 19)
abline(a = 0, b = 1, col = "red", lwd = 2)
Standardized insect counts compared with standard normal quantiles; the red reference is y = x.
The points bend away from the red line, especially at the lower end. The counts do not look normally distributed.
T3
Run a Shapiro-Wilks test on count to check for normality. (0.5 pt) State your hypotheses and elaborate what you deduce from the result. (0.5 pt)
Hypotheses:
H0: the insect counts follow a normal distribution. H1: they do not.
count_normality <- shapiro.test(Sprays$count)
count_normality##
## Shapiro-Wilk normality test
##
## data: Sprays$count
## W = 0.922, p-value = 0.00025The p-value is 0.0002525, which is below 0.05. I reject H0. The combined counts are not normally distributed, which agrees with the plot.
T4
We want to check if the variance in insect counts is the same with the 6 different spray types used. First, using ggplot2 make a boxplot showing insect count (dependent variable) versus all 6 different insecticide sprays (independent variable) (0.5 pts). Make the boxplots filled with green color and set the transparency to 0.3 (0.5 pts). Label your axes appropriately and add a title. Use the black and white theme. (0.5 pts)
ggplot(Sprays, aes(x = spray, y = count)) +
geom_boxplot(fill = "green", alpha = 0.3) +
labs(title = "Insect counts by insecticide spray",
x = "Insecticide spray", y = "Insect count") +
theme_bw()
Insect counts for each of the six insecticide sprays.
Sprays C and E have less spread, while B and F have more spread.
T5
Now run a statistical test to check if the variance in insect counts is the same with the 6 insecticide types (1 pt). State your hypotheses first (0.25 pt), specify why you chose this particular test (0.25 pts) and then discuss the result (0.5 pt)
Hypotheses:
H0: all six sprays have the same population variance. H1: at least one variance is different.
I use the Fligner-Killeen test because it can compare several groups without requiring normal data. I assume the experimental units are independent.
aggregate(count ~ spray, data = Sprays, FUN = var)spray_variance <- fligner.test(count ~ spray, data = Sprays)
spray_variance##
## Fligner-Killeen test of homogeneity of variances
##
## data: count by spray
## Fligner-Killeen:med chi-squared = 14.5, df = 5, p-value = 0.013The p-value is 0.01282, so I reject H0. The variances are not all equal.
T6
Let’s move to the built-in PlantGrowth dataset (Results from an experiment to compare dried weight of plants). Since this is a relatively small dataset (n=30), use Lilliefors Test to check its normality (0.25 pts). Thereafter, based on your results, find out what is the probablity of the plant weight being between 4.5 and 5.5 g. (0.5 pts)
Hypotheses:
H0: plant weights follow a normal distribution. H1: they do not.
I use the Lilliefors test to check normality.
data("PlantGrowth")
plant_normality <- lillie.test(PlantGrowth$weight)
plant_normality##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: PlantGrowth$weight
## D = 0.0934, p-value = 0.72The p-value is 0.7242, which is above 0.05. I fail to reject H0. A normal model seems reasonable, but this does not prove normality.
I use the sample mean and standard deviation to estimate the probability:
plant_mean <- mean(PlantGrowth$weight)
plant_sd <- sd(PlantGrowth$weight)
plant_probability <- pnorm(5.5, mean = plant_mean, sd = plant_sd) -
pnorm(4.5, mean = plant_mean, sd = plant_sd)
data.frame(mean_g = plant_mean, sd_g = plant_sd,
probability = plant_probability)The mean is 5.073 g and the standard deviation is 0.701 g. Using a normal model for all plants together, the chance of a weight between 4.5 and 5.5 g is about 52.18%.
T7
Download from moodle the dataset called “Dineout.csv” and import it into R. This is a made-up dataset of 100 random people asked how much on average they spend when dining out in either Yerevan or Gyumri.
Show the first 8 rows of the dataset below (0.25 pt)
Dineout <- read.csv("dineout.csv", row.names = 1)
Dineout$city <- factor(Dineout$city)
stopifnot(!anyNA(Dineout), is.numeric(Dineout$spending))
head(Dineout, 8)T8
Use the ggplot2 package to make a boxplot with the dependent variable on the y-axis and the independent categorical variable on the x-axis. Label your chart and axes appropriately. (0.75 pt)
ggplot(Dineout, aes(x = city, y = spending)) +
geom_boxplot(fill = "lightblue") +
labs(title = "Dining-out spending in Yerevan and Gyumri",
x = "City", y = "Spending per dining occasion (AMD)") +
theme_bw()
Average dining-out spending reported by respondents in each city.
T9
Test if the average dining cost in both cities is significantly different from 3500 AMD. (2 pts)
Grading details for the above question:
- Choosing and running the appropriate test (0.75 pt)
- Making sure the assumption(s) for your test are met (0.50 pt)
- Stating all the hypotheses correctly (0.5)
- Analyzing and discussing your results correctly (0.25)I use a two-sided one-sample t-test for each city to compare its mean with 3,500 AMD.
For each city: H0: mean spending = 3,500 AMD. H1: mean spending is different from 3,500 AMD.
Spending is numerical, and I assume the randomly selected people are independent. I check normality and outliers for each city. These tests do not require equal variances between cities.
city_values <- split(Dineout$spending, Dineout$city)
city_checks <- do.call(rbind, lapply(names(city_values), function(city) {
x <- city_values[[city]]
data.frame(city = city, n = length(x), mean_AMD = mean(x),
sd_AMD = sd(x), shapiro_p = shapiro.test(x)$p.value,
boxplot_flags = length(boxplot.stats(x)$out),
max_absolute_z = max(abs((x - mean(x)) / sd(x))))
}))
knitr::kable(city_checks, digits = 4, row.names = FALSE)| city | n | mean_AMD | sd_AMD | shapiro_p | boxplot_flags | max_absolute_z |
|---|---|---|---|---|---|---|
| Gyumri | 50 | 4011.8 | 471.35 | 0.8419 | 2 | 2.8238 |
| Yerevan | 50 | 4013.2 | 503.96 | 0.3210 | 0 | 2.5982 |
ggplot(Dineout, aes(sample = spending)) +
stat_qq() + stat_qq_line(colour = "red") +
facet_wrap(~ city) +
labs(title = "Normality checks for dining spending",
x = "Theoretical normal quantiles", y = "Spending (AMD)") +
theme_bw()
Normal Q-Q plots of spending within each city.
Both normality p-values are above 0.05, and the Q-Q plots look roughly normal. Each city has 50 people. Gyumri has two mild outliers, but neither is more than three standard deviations from the mean. The t-tests seem suitable.
city_tests <- lapply(city_values, function(x) {
t.test(x, mu = 3500, alternative = "two.sided", conf.level = 0.95)
})
city_tests## $Gyumri
##
## One Sample t-test
##
## data: x
## t = 7.68, df = 49, p-value = 6e-10
## alternative hypothesis: true mean is not equal to 3500
## 95 percent confidence interval:
## 3877.8 4145.7
## sample estimates:
## mean of x
## 4011.8
##
##
## $Yerevan
##
## One Sample t-test
##
## data: x
## t = 7.2, df = 49, p-value = 3.3e-09
## alternative hypothesis: true mean is not equal to 3500
## 95 percent confidence interval:
## 3869.9 4156.4
## sample estimates:
## mean of x
## 4013.2city_results <- do.call(rbind, lapply(names(city_tests), function(city) {
z <- city_tests[[city]]
data.frame(city = city, mean_AMD = unname(z$estimate),
difference_from_3500 = unname(z$estimate) - 3500,
t = unname(z$statistic), df = unname(z$parameter),
p = z$p.value, CI_low = z$conf.int[1], CI_high = z$conf.int[2])
}))
city_results$Holm_adjusted_p <- p.adjust(city_results$p, method = "holm")
knitr::kable(city_results, digits = c(0, 2, 2, 3, 0, 12, 2, 2, 12),
row.names = FALSE)| city | mean_AMD | difference_from_3500 | t | df | p | CI_low | CI_high | Holm_adjusted_p |
|---|---|---|---|---|---|---|---|---|
| Gyumri | 4011.8 | 511.77 | 7.678 | 49 | 5.960e-10 | 3877.8 | 4145.7 | 1.193e-09 |
| Yerevan | 4013.2 | 513.16 | 7.200 | 49 | 3.251e-09 | 3869.9 | 4156.4 | 3.251e-09 |
The mean is 4013.16 AMD in Yerevan and 4011.77 AMD in Gyumri. Both p-values are below 0.05, so I reject both null hypotheses. Average spending is higher than 3,500 AMD in both cities. Adjusting for the two tests does not change this conclusion.
T10
Download the “fatalities.csv” dataset from moodle and import it here. The data is taken from the World Health Organization (WHO) and reports fatalities (deaths) from car accidents in different countries and years.
Firstly, import the said dataset into R and call it Deaths. Produce the first 7 rows of the dataset. (0.25 pts)
Deaths <- read.csv("fatalities.csv", fileEncoding = "UTF-8-BOM")
stopifnot(!anyNA(Deaths[c("Country", "Year", "Gender", "Deaths")]))
head(Deaths, 7)The dataset includes 7,320 incidents from all countries. For the next step, I will help you pick out from the dataset only the incidents that have happened in Armenia.
To do that, we will use a function called subset to create a new object in R that we’ll call “Deaths.Arm”. It will select from the Deaths dataset you imported only the rows that specify Armenia within the “Country” column:
Deaths.Arm = subset(Deaths, Country == "Armenia")
print(Deaths.Arm)## Indicator ValueType
## 25 Estimated road traffic death rate (per 100 000 population) numeric
## 207 Estimated road traffic death rate (per 100 000 population) numeric
## 577 Estimated road traffic death rate (per 100 000 population) numeric
## 715 Estimated road traffic death rate (per 100 000 population) numeric
## 895 Estimated road traffic death rate (per 100 000 population) numeric
## 1082 Estimated road traffic death rate (per 100 000 population) numeric
## 1244 Estimated road traffic death rate (per 100 000 population) numeric
## 1457 Estimated road traffic death rate (per 100 000 population) numeric
## 1632 Estimated road traffic death rate (per 100 000 population) numeric
## 1799 Estimated road traffic death rate (per 100 000 population) numeric
## 2000 Estimated road traffic death rate (per 100 000 population) numeric
## 2172 Estimated road traffic death rate (per 100 000 population) numeric
## 2369 Estimated road traffic death rate (per 100 000 population) numeric
## 2527 Estimated road traffic death rate (per 100 000 population) numeric
## 2745 Estimated road traffic death rate (per 100 000 population) numeric
## 2869 Estimated road traffic death rate (per 100 000 population) numeric
## 3106 Estimated road traffic death rate (per 100 000 population) numeric
## 3272 Estimated road traffic death rate (per 100 000 population) numeric
## 3465 Estimated road traffic death rate (per 100 000 population) numeric
## 3655 Estimated road traffic death rate (per 100 000 population) numeric
## 3831 Estimated road traffic death rate (per 100 000 population) numeric
## 3999 Estimated road traffic death rate (per 100 000 population) numeric
## 4209 Estimated road traffic death rate (per 100 000 population) numeric
## 4332 Estimated road traffic death rate (per 100 000 population) numeric
## 4409 Estimated road traffic death rate (per 100 000 population) numeric
## 4544 Estimated road traffic death rate (per 100 000 population) numeric
## 4766 Estimated road traffic death rate (per 100 000 population) numeric
## 4921 Estimated road traffic death rate (per 100 000 population) numeric
## 5306 Estimated road traffic death rate (per 100 000 population) numeric
## 5483 Estimated road traffic death rate (per 100 000 population) numeric
## 5663 Estimated road traffic death rate (per 100 000 population) numeric
## 5850 Estimated road traffic death rate (per 100 000 population) numeric
## 6042 Estimated road traffic death rate (per 100 000 population) numeric
## 6213 Estimated road traffic death rate (per 100 000 population) numeric
## 6243 Estimated road traffic death rate (per 100 000 population) numeric
## 6405 Estimated road traffic death rate (per 100 000 population) numeric
## 6606 Estimated road traffic death rate (per 100 000 population) numeric
## 6778 Estimated road traffic death rate (per 100 000 population) numeric
## 7165 Estimated road traffic death rate (per 100 000 population) numeric
## 7306 Estimated road traffic death rate (per 100 000 population) numeric
## Region.Abr Region Country.Abr Country Year Gender Deaths
## 25 EUR Europe ARM Armenia 2019 Female 10.66
## 207 EUR Europe ARM Armenia 2019 Male 30.41
## 577 EUR Europe ARM Armenia 2018 Male 33.62
## 715 EUR Europe ARM Armenia 2018 Female 8.64
## 895 EUR Europe ARM Armenia 2017 Male 25.25
## 1082 EUR Europe ARM Armenia 2017 Female 8.53
## 1244 EUR Europe ARM Armenia 2016 Male 22.94
## 1457 EUR Europe ARM Armenia 2016 Female 9.27
## 1632 EUR Europe ARM Armenia 2015 Male 27.32
## 1799 EUR Europe ARM Armenia 2015 Female 7.37
## 2000 EUR Europe ARM Armenia 2014 Male 26.74
## 2172 EUR Europe ARM Armenia 2014 Female 7.72
## 2369 EUR Europe ARM Armenia 2013 Male 27.92
## 2527 EUR Europe ARM Armenia 2013 Female 6.92
## 2745 EUR Europe ARM Armenia 2012 Male 29.62
## 2869 EUR Europe ARM Armenia 2012 Female 5.92
## 3106 EUR Europe ARM Armenia 2011 Male 28.56
## 3272 EUR Europe ARM Armenia 2011 Female 8.20
## 3465 EUR Europe ARM Armenia 2010 Male 27.75
## 3655 EUR Europe ARM Armenia 2010 Female 9.50
## 3831 EUR Europe ARM Armenia 2009 Male 29.13
## 3999 EUR Europe ARM Armenia 2009 Female 7.95
## 4209 EUR Europe ARM Armenia 2008 Male 29.78
## 4332 EUR Europe ARM Armenia 2008 Female 6.36
## 4409 EUR Europe ARM Armenia 2007 Female 11.13
## 4544 EUR Europe ARM Armenia 2007 Male 24.83
## 4766 EUR Europe ARM Armenia 2006 Female 10.11
## 4921 EUR Europe ARM Armenia 2006 Male 26.58
## 5306 EUR Europe ARM Armenia 2005 Male 28.61
## 5483 EUR Europe ARM Armenia 2005 Female 9.18
## 5663 EUR Europe ARM Armenia 2004 Male 28.36
## 5850 EUR Europe ARM Armenia 2004 Female 9.67
## 6042 EUR Europe ARM Armenia 2003 Male 28.90
## 6213 EUR Europe ARM Armenia 2003 Female 9.26
## 6243 EUR Europe ARM Armenia 2002 Female 10.93
## 6405 EUR Europe ARM Armenia 2002 Male 28.29
## 6606 EUR Europe ARM Armenia 2001 Female 10.43
## 6778 EUR Europe ARM Armenia 2001 Male 29.69
## 7165 EUR Europe ARM Armenia 2000 Male 32.16
## 7306 EUR Europe ARM Armenia 2000 Female 8.57unique(Deaths$Indicator)## [1] "Estimated road traffic death rate (per 100 000 population)"The file has 7320 rows. Its values are death rates per 100,000 people, not numbers of deaths. I use those units in the graphs.
T11
Here is our research question: Based on the number of fatalities (Deaths) due to car accidents, we want to test if men in Armenia are indeed better drivers than women as the stereotype goes. Assuming the data follows a normal distribution, state your hypotheses and test them, making sure you satisfy all other assumptions (1.75 pt)
Grading details for the above question:
- Stating all the hypotheses correctly (0.5)
- Assume normality. Make sure the rest of the assumption(s) for your test are met (0.25 pt)
- Choose and run the appropriate test (0.75 pt)
- Analyze and discussing your results correctly (0.25)I test whether men have a lower death rate than women. The difference is the male rate minus the female rate in the same year.
H0: the mean difference is zero or positive. H1: the mean difference is negative, meaning a lower male rate.
I use a one-sided paired t-test because male and female rates are matched by year. I assume normal differences as the question asks. Equal variances are not needed.
stopifnot(all(Deaths.Arm$Gender %in% c("Female", "Male")))
arm_pairs <- with(Deaths.Arm, table(Year, Gender))
arm_pairs## Gender
## Year Female Male
## 2000 1 1
## 2001 1 1
## 2002 1 1
## 2003 1 1
## 2004 1 1
## 2005 1 1
## 2006 1 1
## 2007 1 1
## 2008 1 1
## 2009 1 1
## 2010 1 1
## 2011 1 1
## 2012 1 1
## 2013 1 1
## 2014 1 1
## 2015 1 1
## 2016 1 1
## 2017 1 1
## 2018 1 1
## 2019 1 1stopifnot(all(arm_pairs == 1))
Arm.wide <- reshape(Deaths.Arm[c("Year", "Gender", "Deaths")],
idvar = "Year", timevar = "Gender", direction = "wide")
Arm.wide <- Arm.wide[order(Arm.wide$Year), ]
stopifnot(!anyNA(Arm.wide))
Arm.wide$difference <- Arm.wide$Deaths.Male - Arm.wide$Deaths.Female
data.frame(pairs = nrow(Arm.wide),
mean_male_rate = mean(Arm.wide$Deaths.Male),
mean_female_rate = mean(Arm.wide$Deaths.Female),
mean_difference = mean(Arm.wide$difference))boxplot.stats(Arm.wide$difference)$out## numeric(0)There are 20 complete year pairs and no boxplot outliers in the differences. The test assumes years are independent. Since nearby years may be related, this assumption is a limitation.
arm_test <- t.test(Arm.wide$Deaths.Male, Arm.wide$Deaths.Female,
paired = TRUE, alternative = "less", mu = 0)
arm_test##
## Paired t-test
##
## data: Arm.wide$Deaths.Male and Arm.wide$Deaths.Female
## t = 28.5, df = 19, p-value = 1
## alternative hypothesis: true mean difference is less than 0
## 95 percent confidence interval:
## -Inf 20.689
## sample estimates:
## mean difference
## 19.507The male rate is higher by 19.507 deaths per 100,000 on average. The one-sided p-value is about 1.000, so I fail to reject H0. There is no support for a lower male death rate.
This does not tell us who drives better. Deaths may include passengers and pedestrians, and the data do not show how much people drive or who caused each accident.
T12
Using ggplot2, make a line graph of fatal accidents in Armenia vs years. Make sure the dependent variable is on the y-axis (0.5 pts). Use different colored points to distinguish the 2 genders (0.25 pts) and make the line size =3 (0.25 pts). Make sure you label your axes and title the graph appropriately (0.25 pts)
ggplot(Deaths.Arm, aes(x = Year, y = Deaths, colour = Gender,
group = Gender)) +
geom_line(linewidth = 3) +
geom_point(size = 2.5) +
scale_colour_manual(values = c(Female = "#C23B68", Male = "#2475B0")) +
labs(title = "Road traffic death rates in Armenia",
x = "Year", y = "Deaths per 100,000 population", colour = "Sex") +
theme_bw()
Estimated road traffic death rates in Armenia, by year and sex.
The male death rate is higher in every year shown.
T13
Going back to the original dataset you imported (Deaths), using ggplot2 make a boxplot of deaths per region (0.5 pt). Use the dark theme (0.25) and fill the boxplots with orange (0.25).
ggplot(Deaths, aes(x = Region, y = Deaths)) +
geom_boxplot(fill = "orange") +
labs(title = "Road traffic death rates by region",
x = "Region", y = "Deaths per 100,000 population") +
theme_dark() +
theme(axis.text.x = element_text(angle = 30, hjust = 1))
Distribution of sex-specific country-year road traffic death rates in each region.
Each value is a rate for one sex in one country and year.
T14
Similar to the approach above, make a subset from the Deaths dataset that selects out data from Europe (Region) alone. Call the new object “Deaths.Europe” (0.25 pts).
Deaths.Europe <- subset(Deaths, Region == "Europe")
head(Deaths.Europe)nrow(Deaths.Europe)## [1] 2000T15
Using the new dataset created, check what is the average of fatalities (Deaths) in Europe (0.5 pts). Important to note that the data is split by gender within the dataset, so for the whole population, you’ll need to account for both genders. Following the question, I add the male and female values for each country and year, then find the average.
europe_pairs <- with(Deaths.Europe, table(Country, Year, Gender))
stopifnot(all(Deaths.Europe$Gender %in% c("Female", "Male")),
all(europe_pairs == 1))
Europe.yearly <- aggregate(Deaths ~ Country + Year,
data = Deaths.Europe, FUN = sum)
names(Europe.yearly)[3] <- "sum_of_sex_specific_values"
europe_assignment_mean <- mean(Europe.yearly$sum_of_sex_specific_values)
data.frame(countries = length(unique(Europe.yearly$Country)),
country_years = nrow(Europe.yearly),
assignment_mean = europe_assignment_mean)This calculation gives an average of 22.3309.
However, these values are rates, so adding them does not give a true total death count or overall rate. We would need the male and female population sizes to calculate that correctly.
T16
Accordingly, determine what is the probability of deaths being higher than 60 in any given year and country in Europe (1 point). Hint to help you: This question deals with a frequency of an event (deaths) that occurs within a specific interval (year).
Following the hint, I use a Poisson model with the average from T15. More than 60 means 61 or more.
lambda_assignment <- europe_assignment_mean
probability_over_60 <- ppois(60, lambda = lambda_assignment,
lower.tail = FALSE)
data.frame(lambda = lambda_assignment,
probability_more_than_60 = probability_over_60)The calculated probability is 1.1788e-11, which is almost zero.
This follows the exercise’s count-based method. It is not a reliable real-world probability because the file contains rates, not counts, and countries have different populations.
T17
Import the “treatment.csv” dataset which includes differences in blood sugar concentration after the administration of two different drugs (expressed in difference before and after drug use). To minimize the impact of other factors (e.g. genetic, lifestyle, etc.), each of the participants of this study were tested twice. Once with drug 1 and the second time with drug 2 after a few days. We want to check if the drugs differ from each other in efficacy (i.e. one is significantly more effective than the other). Run the appropriate test, after checking all the assumptions. (1.75 pts)
Grading details for the above question:
- Stating all the hypotheses correctly (0.25)
- Making sure the assumption(s) for your test are met (0.5 pt)
- Choosing and running the appropriate test (0.75 pt)
- Analyzing and discussing your results correctly (0.25)Hint to help you: In the above dataset where you have the 2 groups arranged separately in different columns (and you don’t have another column for categories / factor), for variance and t-testing, do not use the ~ sign (that works only with factors), instead use a simple comma (,) between the two groups.
The same people tried both drugs, so the data are paired. I use drug 1 minus drug 2 for each person. H0: the mean difference is zero. H1: the mean difference is not zero.
Treatment.raw <- read.csv("treatment.csv", check.names = FALSE)
# Remove only entirely empty trailing rows; do not drop genuine measurements.
empty_rows <- rowSums(!is.na(Treatment.raw)) == 0
Treatment <- Treatment.raw[!empty_rows, ]
stopifnot(!anyNA(Treatment), !anyDuplicated(Treatment$ID),
is.numeric(Treatment[["drug 1"]]),
is.numeric(Treatment[["drug 2"]]))
drug1 <- Treatment[["drug 1"]]
drug2 <- Treatment[["drug 2"]]
drug_difference <- drug1 - drug2
data.frame(empty_rows_removed = sum(empty_rows),
complete_pairs = nrow(Treatment),
mean_drug1 = mean(drug1), mean_drug2 = mean(drug2),
mean_difference = mean(drug_difference),
median_difference = median(drug_difference))After removing 19 empty rows, there are 20 complete pairs. I check normality of the differences. Equal variances are not needed. I assume people are independent and the first drug has worn off before the second test; the file cannot confirm this.
drug_normality <- shapiro.test(drug_difference)
drug_normality##
## Shapiro-Wilk normality test
##
## data: drug_difference
## W = 0.803, p-value = 0.00096boxplot.stats(drug_difference)$out## [1] -4.6 -4.6old_par <- par(mfrow = c(1, 2))
qqnorm(drug_difference, main = "Q-Q plot of paired differences", pch = 19)
qqline(drug_difference, col = "red", lwd = 2)
boxplot(drug_difference, main = "Paired differences",
ylab = "Drug 1 change minus drug 2 change", col = "lightblue")
abline(h = 0, lty = 2)
Normality and outlier checks for the paired drug differences.
par(old_par)For normality, H0: the differences are normal; H1: they are not. The p-value is 0.0009573, so I reject normality. There are also two outliers at -4.6. This makes the paired t-test less suitable.
I mainly use an exact paired sign test because it does not need normal or symmetric differences. H0: either drug is equally likely to give the larger change among non-ties. H1: one is more likely. This checks the direction of the difference, not the mean.
# Round to the recorded precision to avoid floating-point comparisons at zero.
d <- round(drug_difference, 10)
nonzero_d <- d[d != 0]
sign_test <- binom.test(sum(nonzero_d > 0), length(nonzero_d),
p = 0.5, alternative = "two.sided")
data.frame(drug1_greater = sum(d > 0), drug2_greater = sum(d < 0),
tied = sum(d == 0))sign_test##
## Exact binomial test
##
## data: sum(nonzero_d > 0) and length(nonzero_d)
## number of successes = 0, number of trials = 18, p-value = 7.6e-06
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
## 0.0000 0.1853
## sample estimates:
## probability of success
## 0Drug 2 gives a larger change for all 18 non-tied people; 2 people have equal changes. The p-value is 7.629e-06, so I reject H0. Drug 2’s recorded change is larger by 1.58 units on average.
I also show the paired t-test and Wilcoxon test for comparison:
drug_t_test <- t.test(drug1, drug2, paired = TRUE,
alternative = "two.sided")
drug_t_test##
## Paired t-test
##
## data: drug1 and drug2
## t = -5.9, df = 19, p-value = 1.1e-05
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
## -2.1403 -1.0197
## sample estimates:
## mean difference
## -1.58# Differences are rounded before rank testing to preserve recorded ties.
drug_wilcoxon <- wilcox.test(round(drug1, 10), round(drug2, 10),
paired = TRUE, alternative = "two.sided",
exact = FALSE, correct = TRUE)
drug_wilcoxon##
## Wilcoxon signed rank test with continuity correction
##
## data: round(drug1, 10) and round(drug2, 10)
## V = 0, p-value = 0.00021
## alternative hypothesis: true location shift is not equal to 0The paired t-test (p = 1.106e-05) and Wilcoxon test (p = 0.0002104) also find a difference. However, normality fails and the differences look uneven around the center, so I rely more on the sign test.
If a positive change means a blood-sugar reduction, drug 2 works better. Otherwise, I can only say it gives a larger recorded change. The differences repeat in the second set of ten people, so the result assumes these are still separate, independent participants.