Install and load the packages you need to produce the report here:
install.packages("dplyr") #data manipulation
install.packages("ggplot2") #data visualisations
install.packages("tidyverse")
library(dplyr)
library(ggplot2)
library(knitr)
library(tidyverse)
In this section of the assessment, there are two groups, referred to as group A and group B, Group A receives the flu shot, and group B receives a placebo shot. The contagion rate for group A (those with the jab) is 10%, and for group B it was 30%.
#Setting seed for reproductability
set.seed(3388921)
# Setting parameters for Binomial distribution
n_group_a <- 100
n_group_b <- 100
p_a <- 0.1
p_b <- 0.3
# Generate synthetic data
flu_cases_group_a <- rbinom(1, n_group_a, p_a)
flu_cases_group_b <- rbinom(1, n_group_b, p_b)
# Calculate the number of people who did not contract the flu
no_flu_cases_group_a <- n_group_a - flu_cases_group_a
no_flu_cases_group_b <- n_group_b - flu_cases_group_b
# Making a data frame
df <- data.frame(
Condition = c("Contracted the flu", "Did not contract the flu"),
Group_A = c(flu_cases_group_a, no_flu_cases_group_a),
Group_B = c(flu_cases_group_b, no_flu_cases_group_b)
)
print(df)
## Condition Group_A Group_B
## 1 Contracted the flu 10 30
## 2 Did not contract the flu 90 70
| Condition Group | Group A | Group B |
|---|---|---|
| Contracted the flu | 10 | 30 |
| Did not contract the flu | 90 | 70 |
| Total: | 100 | 100 |
What is your estimated probability that a person receiving the placebo will not contract the flu during the winter season?
When looking to the table, we can see that the total number of people who received the placebo jab was 100, with 70 of them not contracting the flu, thus leaving the estimated probability at 0.7. However with the additional caveat of this being during the winter season, we would need to simply multiply 70/100 by the relevant variable (k) (in this instance there is no specified impact towards the participants and thus I will leave it as no impact). In conclusion for this section the estimated probability remains at 0.7
Derive a 95% confidence interval for the proportion of people receiving the placebo who do not contract the flu in that winter season. Give a non-technical explanation of your result.
install.packages("binom")
library(binom)
#calculating total number of people in group B.
total_group_b <- sum(df$Group_B)
cat("Total number of people in group B:", total_group_b, "\n")
## Total number of people in group B: 100
#calculating total number of people in group B who did not contract the flu
p_hat <- no_flu_cases_group_b / n_group_b
#using binom confit function to find confidence value.
confint <- binom.confint(no_flu_cases_group_b, n_group_b, conf.level = 0.95, methods = "exact")
# Printing the 95% confidence interval
cat("people in the placebo group who did not contract the flu:", confint$mean, "\n")
## people in the placebo group who did not contract the flu: 0.7
cat("95% Confidence Interval for the placebo group:", confint$lower, "to", confint$upper, "\n")
## 95% Confidence Interval for the placebo group: 0.6001853 to 0.7875936
I underwent conducting a study where 100 people received a placebo flu shot instead of the actual flu vaccine. Approximately 70% of the people who received the placebo did not catch the flu during the winter season. To understand how precise this estimate of 70% is, I formulated a 95% confidence interval. This interval helps me determine that i can be 95% sure or confident that the true proportion of people who would not catch the flu, if everyone in the population received the placebo, lies between approximately 61% and 79%. This range gives us an idea of the variability and reliability of our estimate.
People in the placebo group who did not contract the flu: 0.7
95% Confidence Interval for the placebo group: 0.6001853 to 0.7875936
If 40% of the wider adult population receive the flu vaccine, what percentage would you anticipate contracting the flu that year?
40% of the population are given the flu vaccine, therefore 40% of the total population have received the flu vaccine and are apart of group A and the remainder are given the placebo jab and are part of group B. The contraction rate for those vaccinated is 0.1(10%) and therefor the percentage rate for contracting the flu would be 22%.
# probabilities of contracting the flu if/if not vaccinated.
p_vaccine <- 0.4 # Proportion of population that IS vaccinated
p_no_vaccine <- 0.6 # Proportion of population that IS NOT vaccinated
p_flu_vaccine <- 0.1 # Probability of contracting flu IF vaccinated
p_flu_no_vaccine <- 0.3 # Probability of contracting flu IF NOT vaccinated
# Calculating the overall probability of contracting the flu.
p_flu <- (p_flu_vaccine * p_vaccine) + (p_flu_no_vaccine * p_no_vaccine)
# Converting the rate to percentage.
p_flu_percentage <- p_flu * 100
cat("Anticipated percentage of the wider adult population contracting the flu:", p_flu_percentage, "%\n")
## Anticipated percentage of the wider adult population contracting the flu: 22 %
"Anticipated percentage of the wider adult population contracting the flu: 22 %"
What is your estimate of the probability that a person who is vaccinated against the flu for ten years running, will contract the flu on at least three of those years?
In the second line of code, I am calculating the cumulative probability of contracting the flu for 0, 1, or 2 times out of 10 years. In the third line it will give the probability of contracting the flu at least 3 times out of 10 years. The result we get from this is an estimated probability of 0.0702 (or 7.02%)
# setting the probability and number of years
p_flu_vaccine <- 0.1
n_years <- 10
p_less_than_3 <- pbinom(2, n_years, p_flu_vaccine)
p_at_least_3 <- 1 - p_less_than_3
cat("Estimated probability of contracting the flu on at least three of the ten years:", p_at_least_3, "\n")
## Estimated probability of contracting the flu on at least three of the ten years: 0.07019083
This question explores the distribution of the sample mean when the underlying variable has an exponential distribution with mean 10.
# Defining the mean of the exponential distribution
mean_exp <- 10
# Calculating the rate
rate <- 1 / mean_exp
# Producing 1000 observations from the exponential distribution
observations <- rexp(1000, rate)
# Display the first few observations
head(observations)
## [1] 12.288405 1.867218 11.856425 11.765654 5.871312 13.055570
[1] 12.288405 1.867218 11.856425 11.765654 5.871312 13.055570
Compare the sample estimates for the mean and standard deviation with the population mean and standard deviation for the exponential distribution. Discuss.
#Producing the smaple mean and sd.
sample_mean <- mean(observations)
sample_sd <- sd(observations)
#Printing sample mean and sd.
cat("Sample mean:", sample_mean, "\n")
## Sample mean: 8.954161
cat("Sample standard deviation:", sample_sd, "\n")
## Sample standard deviation: 8.989938
#Setting population parameters.
population_mean <- 1 / rate
population_sd <- 1 / rate
#Printing the population mean and standard deviation.
cat("Population mean:", population_mean, "\n")
## Population mean: 10
cat("Population standard deviation:", population_sd, "\n")
## Population standard deviation: 10
Sample mean: 8.954161
Population mean: 10
Sample sd: 8.989938
Population sd: 10
Compare the histogram of the sample data with the density function of the exponential distribution. Discuss.
The results from the sample mean and sd match closely to the population mean and sd, which can be reflected in the visualisation above. You see the the bars closely conform to, as well as match the red curve on the plot, and this is because the sample accurately represents the exponential distribution.
Display a histogram of the sample means and discuss its characteristics. Generate 1000 sample means of sample size 2 where the observations are drawn at random from the exponential distribution.
n_samples <- 1000
sample_size <- 2
sample_means <- replicate(n_samples, mean(rexp(sample_size, rate)))
head(sample_means)
## [1] 20.837843 8.561333 22.669389 35.925186 10.408214 8.843121
Showing first few samples.
[1] 20.837843 8.561333 22.669389 35.925186 10.408214 8.843121
Display a histogram of the sample means and discuss its characteristics.
Generate 1000 sample means of sample size 30 where the observations are drawn at random from the exponential distribution.
n_samples <- 1000
sample_size <- 30
sample_means <- replicate(n_samples, mean(rexp(sample_size, rate)))
head(sample_means)
## [1] 8.368059 10.019416 10.786604 12.302778 9.193607 11.230635
Showing first few samples.
[1] 8.368059 10.019416 10.786604 12.302778 9.193607 11.230635
In preparation for your meeting with management, conduct a one-tailed test of the difference between the means, using a 5% significance level. What can you conclude from the test? [Hint: data is normally distributed].
mean_current <- 7.5
mean_new <- 8.2
sd_diff <- 1.9
n <- 20
# Mean difference
mean_diff <- mean_new - mean_current
# Test statistic
t_stat <- mean_diff / (sd_diff / sqrt(n))
# Degrees of freedom
dof <- n - 1
# P-value for one tailed test.
p_value <- pt(t_stat, dof, lower.tail = FALSE)
# Critical value for one tailed test at 5% significance level
alpha <- 0.05
t_critical <- qt(alpha, dof, lower.tail = FALSE)
cat("Test Statistic (t):", t_stat, "\n")
## Test Statistic (t): 1.647629
cat("Critical Value (t_critical):", t_critical, "\n")
## Critical Value (t_critical): 1.729133
cat("P-value:", p_value, "\n")
## P-value: 0.05793374
# Decision
if (t_stat > t_critical) {
cat("Reject the null hypothesis: The new design is significantly more attractive.\n")
} else {
cat("Fail to reject the null hypothesis: There is no significant difference in attractiveness.\n")
}
## Fail to reject the null hypothesis: There is no significant difference in attractiveness.
Test Statistic (t): 1.647629
Critical Value (t_critical): 1.729133
P-value: 0.05793374
Mean difference: 0.7
Looking at the means of the new and old designs, we can see initially that both means are roughly similar in value. This was acheived through the one tailed test that I utilised in the above code. The aim of the one tailed test allows for us to impose a higher level of significance via the critical value onto one object over the other (in this case it was towards the newer design). The results show that despite the that there is little difference between the old and new design. This indicates that the opinion of the new and old design are roughly the same. The target of the new design is to have a greater attractiveness (in this case, a signficantly higher mean) than the older design.
cat("Test Statistic (t):", t_stat, "\n")
## Test Statistic (t): 1.647629
cat("Critical Value (t_critical):", t_critical, "\n")
## Critical Value (t_critical): 1.729133
cat("P-value:", p_value, "\n")
## P-value: 0.05793374
# Decision
if (t_stat > t_critical) {
cat("The new design is significantly more attractive.\n")
} else {
cat("There is no significant difference in attractiveness.\n")
}
## There is no significant difference in attractiveness.
"There is no significant difference in attractiveness."
Advise to management:
Currently, my advise to the company/management would be to not
proceed with this new design. I will explain my reasoning as well as
produce some alternative margins that the company can aim towards, as
this will assist with knowing when it becomes feesable and worthy to
make the change over to a new design. As you can see above, the output
that was produced deemed that “there is no significant difference in
attractiveness”. Thus when factoring in costs and all other re-branding
initiatives required to make the switch to a newer design, it simply
isnt worth it, as the design isn’t much different from the older
one.
Even with the lack of difference in attractiveness, it still will cost
the company the same amount as if it was significantly different.
Therefore it is better for the company to have a design that is
marginally better and pay for the re-branding then. We can produce a
margin which can allow for the company to know when it is worth while
considering the change over.
alpha <- 0.08
t_critical <- qt(alpha, dof, lower.tail = FALSE)
# Output the results
cat("Test Statistic (t):", t_stat, "\n")
## Test Statistic (t): 1.647629
cat("Critical Value (t_critical):", t_critical, "\n")
## Critical Value (t_critical): 1.462314
cat("P-value:", p_value, "\n")
## P-value: 0.05793374
# Decision
if (t_stat > t_critical) {
cat("The new design is significantly more attractive.\n")
} else {
cat("There is no significant difference in attractiveness.\n")
}
## The new design is significantly more attractive.
With the P value being changed from 0.05 to 0.08 we get the output "The new design is significantly more attractive."
Generate 10 pairs of observations from a bivariate Normal distribution with parameter values
(𝜇1, 𝜇2, 𝜎1, 𝜎2, 𝜌) = (50, 55, 10, 10, 0.8).
set.seed(3388921)
library(MASS) #to be able to utilise mvrnorm function.
mu1 <- 50
mu2 <- 55
sigma1 <- 10
sigma2 <- 10
rho <- 0.8
mean <- c(mu1, mu2)
cov <- matrix(c(sigma1^2, rho * sigma1 * sigma2, rho * sigma1 * sigma2, sigma2^2), nrow = 2)
observations <- MASS::mvrnorm(n = 10, mu = mean, Sigma = cov)
observations
## [,1] [,2]
## [1,] 50.68642 55.21542
## [2,] 61.33849 65.48350
## [3,] 51.77141 57.71016
## [4,] 55.33727 60.19133
## [5,] 47.57414 51.17231
## [6,] 47.20112 53.14417
## [7,] 38.91472 47.94533
## [8,] 43.15137 48.64500
## [9,] 52.65786 58.09026
## [10,] 64.56410 70.29058
[,1] [,2]
[1,] 50.68642 55.21542
[2,] 61.33849 65.48350
[3,] 51.77141 57.71016
[4,] 55.33727 60.19133
[5,] 47.57414 51.17231
[6,] 47.20112 53.14417
[7,] 38.91472 47.94533
[8,] 43.15137 48.64500
[9,] 52.65786 58.09026
[10,] 64.56410 70.29058
Decide on an appropriate hypothesis test of the difference between the two means based only on these 10 observations (i.e. as though you are unaware of the parameter values). Explain how you arrived at that decision.
I am working with only ten observations, so the outlook is limited. Because of this limitation, there could be a potential problem of not being able to detect any errors. But becuase I am assuming the observations were paired appropriately, I have decided to go with the paired t-test method regardless. The aim is to determine if there is a significant difference between the two means (from the 10 observations).
Conduct the test using a 5% significance level and write up your findings. Conclude this part of the question by reflecting on the ability of the test to reach a correct conclusion.
observations <- matrix(c(53.64, 62.00, 63.21, 52.33, 58.48, 38.43, 73.61, 39.28, 61.49, 64.23,
54.58, 70.65, 64.38, 59.47, 59.48, 47.76, 82.34, 56.73, 61.74, 70.13), ncol = 2, byrow = TRUE)
# paired t-test
result <- t.test(observations[,1], observations[,2], paired = TRUE)
result
##
## Paired t-test
##
## data: observations[, 1] and observations[, 2]
## t = 1.3917, df = 9, p-value = 0.1974
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
## -4.499313 18.887313
## sample estimates:
## mean difference
## 7.194
Paired t-test
data: observations[, 1] and observations[, 2]
t = 1.3917, df = 9, p-value = 0.1974
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-4.499313 18.887313
sample estimates:
mean of the differences
7.194
P value is 0.1974.
The p value is greater than 0.5. which indicates that there is no significant difference.
Generate 30 pairs of observations from a bivariate Normal distribution with the same parameter values, (𝜇1, 𝜇2, 𝜎1, 𝜎2, 𝜌) = (50, 55, 10, 10, 0.8).
set.seed(3388921)
library(MASS)
mu1 <- 50
mu2 <- 55
sigma1 <- 10
sigma2 <- 10
rho <- 0.8
mean <- c(mu1, mu2)
covariance_matrix <- matrix(c(sigma1^2, rho * sigma1 * sigma2,
rho * sigma1 * sigma2, sigma2^2), nrow = 2)
observations <- mvrnorm(n = 30, mu = mean, Sigma = covariance_matrix)
print(observations)
## [,1] [,2]
## [1,] 51.97208 53.92976
## [2,] 55.41466 71.40733
## [3,] 54.03341 55.44815
## [4,] 59.75778 55.77083
## [5,] 48.43334 50.31311
## [6,] 49.45759 50.88770
## [7,] 47.55109 39.30895
## [8,] 47.21687 44.57950
## [9,] 54.13162 56.61650
## [10,] 63.78845 71.06623
## [11,] 54.75683 48.83017
## [12,] 52.08862 50.34639
## [13,] 51.99022 55.82604
## [14,] 52.59800 51.96418
## [15,] 50.00769 50.78683
## [16,] 54.29155 53.53760
## [17,] 58.95529 58.13654
## [18,] 46.16423 60.31668
## [19,] 52.92841 53.36878
## [20,] 52.17890 55.00056
## [21,] 52.59980 57.72016
## [22,] 46.43819 46.84452
## [23,] 50.61884 47.81712
## [24,] 39.48144 54.12659
## [25,] 68.01743 66.06724
## [26,] 31.36079 42.35799
## [27,] 45.69128 48.87731
## [28,] 39.30496 47.90610
## [29,] 51.20918 63.55068
## [30,] 45.45120 42.63101
[,1] [,2]
[1,] 51.97208 53.92976
[2,] 55.41466 71.40733
[3,] 54.03341 55.44815
[4,] 59.75778 55.77083
[5,] 48.43334 50.31311
[6,] 49.45759 50.88770
[7,] 47.55109 39.30895
[8,] 47.21687 44.57950
[9,] 54.13162 56.61650
[10,] 63.78845 71.06623
[11,] 54.75683 48.83017
[12,] 52.08862 50.34639
[13,] 51.99022 55.82604
[14,] 52.59800 51.96418
[15,] 50.00769 50.78683
[16,] 54.29155 53.53760
[17,] 58.95529 58.13654
[18,] 46.16423 60.31668
[19,] 52.92841 53.36878
[20,] 52.17890 55.00056
[21,] 52.59980 57.72016
[22,] 46.43819 46.84452
[23,] 50.61884 47.81712
[24,] 39.48144 54.12659
[25,] 68.01743 66.06724
[26,] 31.36079 42.35799
[27,] 45.69128 48.87731
[28,] 39.30496 47.90610
[29,] 51.20918 63.55068
[30,] 45.45120 42.63101
Decide on an appropriate hypothesis test of the difference between the two means based only on these 30 observations (i.e. as though you are unaware of the parameter values). Explain how you arrived at that decision.
In this instance I have a larger pool of observations to work with (30) and so with that comes a better level of assurrance of the testing method I choose. Again I decided to go with the paired t-test method as its two sets we are working with, just like last time except that the only difference here is the number of observations has increased. For this reason i've decided to use the same test again.
#t_test_result <- t.test(x, y, paired = TRUE)
#print(t_test_result)
Paired t-test
data: x and y
t = -6.3546, df = 29, p-value = 6.049e-07
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-8.100172 -4.155649
sample estimates:
mean of the differences
-6.127911
Given the p value on the second round of observations (with 30) was below 0.05, meaning there was a significant difference between the two means, in comparison to the first round of observations (with 10) which had a p-value above 0.05, we can conclude that there was strong evidence from the set of 30 observations against the null hypothesis. I believe because this set is more robust (given its larger number of observations) that this would be more correct than from part A’s observations.
Write your plain text here.
set.seed(3388921)
mean_scenario1 <- 50
mean_scenario2 <- 55
mean_scenario3 <- 60
standard_deviation <- 30
num_stores <- 6
num_months <- 6
sales_scenario1 <- matrix(rnorm(num_stores * num_months, mean = mean_scenario1, sd = standard_deviation), nrow = num_months)
sales_scenario2 <- matrix(rnorm(num_stores * num_months, mean = mean_scenario2, sd = standard_deviation), nrow = num_months)
sales_scenario3 <- matrix(rnorm(num_stores * num_months, mean = mean_scenario3, sd = standard_deviation), nrow = num_months)
cat("Scenario 1:\n")
## Scenario 1:
print(sales_scenario1)
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 51.42594 21.31821 54.45289 52.05103 95.9869003 35.568982
## [2,] 84.50359 29.12324 49.30774 53.44603 0.5400581 102.142817
## [3,] 57.08597 59.08857 43.35055 58.41159 33.5065005 32.993640
## [4,] 66.64718 97.20440 54.47328 31.47334 21.8732157 7.371144
## [5,] 40.11227 47.76586 69.11887 39.62134 65.4316966 35.199457
## [6,] 42.64026 45.94439 52.34151 31.98771 23.2506149 33.066562
cat("\n")
cat("Scenario 2:\n")
## Scenario 2:
print(sales_scenario2)
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] -7.812994 49.47781 33.37180 22.03237 49.53241 55.16582
## [2,] 18.772752 28.27645 44.66722 83.44720 42.28713 73.62802
## [3,] 43.069748 34.97868 55.57090 46.39558 46.26992 38.96244
## [4,] 65.804494 27.70665 33.21033 72.08172 45.31336 29.79813
## [5,] 3.170311 27.39925 17.99316 89.82380 30.78141 43.53542
## [6,] 23.018774 98.41388 100.75099 17.90559 76.84616 87.84977
cat("\n")
cat("Scenario 3:\n")
## Scenario 3:
print(sales_scenario3)
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 39.39623 115.88295 78.813863 73.59204 74.26595 118.43552
## [2,] 43.61900 30.74479 106.203392 96.74914 100.46837 53.04286
## [3,] -11.73572 79.33665 57.435971 11.78110 60.14043 122.81086
## [4,] 105.19876 54.46349 5.372779 19.56317 41.57027 38.97824
## [5,] 30.26463 76.14615 96.643433 100.38568 34.78955 35.18297
## [6,] 93.74710 61.07491 -22.076199 21.71300 72.65503 94.17093
Scenario 1:
Month 1: [49, 44, 61, 67, 56, 57]
Month 2: [62, 52, 44, 45, 50, 71]
Month 3: [38, 59, 53, 69, 42, 56]
Month 4: [54, 64, 50, 44, 45, 56]
Month 5: [55, 48, 38, 62, 48, 64]
Month 6: [49, 50, 46, 53, 59, 55]
Scenario 2:
Month 1: [57, 59, 63, 51, 57, 63]
Month 2: [63, 67, 48, 59, 64, 48]
Month 3: [53, 58, 51, 66, 66, 54]
Month 4: [65, 53, 61, 56, 53, 56]
Month 5: [64, 56, 66, 56, 53, 63]
Month 6: [65, 53, 52, 62, 53, 55]
Scenario 3:
Month 1: [64, 63, 74, 61, 64, 69]
Month 2: [64, 60, 63, 68, 75, 57]
Month 3: [52, 65, 72, 71, 59, 66]
Month 4: [55, 67, 70, 71, 71, 67]
Month 5: [59, 69, 66, 67, 63, 64]
Month 6: [69, 72, 72, 68, 58, 71]
Decide on an appropriate test of the difference between the means based only on the observations (i.e. as though you are unaware of the parameter values). Explain how you arrived at that decision.
There are three different scenarios to compare on the sales of each scenario. Because we specifically have three i have decided to go with the one way ANOVA method to compare the difference between means of the three scenarios, as ANOVA is specifically designed for comparing three or more groups at a time and thus will help me see if theres a significant difference in sales between the scenarios.The aim is to see whether or not theres a difference in sales, and if so, what could be contributing to this?
set.seed(3388921)
library(car)
num_scenarios <- 3
num_stores_per_scenario <- 6
num_months <- 6
means <- c(50, 55, 60)
std_dev <- 30
sales_data <- list()
for (i in 1:num_scenarios) {
scenario_sales <- matrix(rnorm(num_stores_per_scenario * num_months, mean = means[i], sd = std_dev), nrow = num_stores_per_scenario, ncol = num_months)
sales_data[[paste("Scenario", i)]] <- scenario_sales
}
data <- data.frame(
Sales = as.vector(unlist(sales_data)),
Scenario = factor(rep(paste("Scenario", 1:num_scenarios), each = num_stores_per_scenario * num_months))
)
by(data$Sales, data$Scenario, shapiro.test)
## data$Scenario: Scenario 1
##
## Shapiro-Wilk normality test
##
## data: dd[x, ]
## W = 0.95805, p-value = 0.1871
##
## ------------------------------------------------------------
## data$Scenario: Scenario 2
##
## Shapiro-Wilk normality test
##
## data: dd[x, ]
## W = 0.96123, p-value = 0.2345
##
## ------------------------------------------------------------
## data$Scenario: Scenario 3
##
## Shapiro-Wilk normality test
##
## data: dd[x, ]
## W = 0.97137, p-value = 0.464
leveneTest(Sales ~ Scenario, data = data)
Conduct the test using a 5% significance level and write up your findings. What can you conclude?
| Df | Sum Sq | Mean Sq | F Value | Pr(>F) | |
|---|---|---|---|---|---|
| Scenario | 2 | 5196 | 2597.8 | 2.974 | 0.0554 |
| Residuals | 105 | 91726 | 873.6 |
The findings show us that the p value is 0.0554 (5.54%). This is just a little higher than the common significance level (5%, we used in the previous question). The p-value here suggests that the difference in mean is insignificant, however it does suggest that it is close to being significant (given its only 0.54% above the threshold). So this does leave us with an indication that promotions can have some effect.
Repeat part (a) under the changed assumption the standard deviation of sales is 25.
set.seed(3388921)
num_scenarios <- 3
num_stores_per_scenario <- 6
num_months <- 6
means <- c(50, 55, 60)
std_dev <- 25 #standard deviation has been changed to 25
sales_data <- list()
for (i in 1:num_scenarios) {
scenario_sales <- matrix(rnorm(num_stores_per_scenario * num_months, mean = means[i], sd = std_dev),
nrow = num_stores_per_scenario, ncol = num_months)
sales_data[[paste("Scenario", i)]] <- scenario_sales
}
data <- data.frame(
Sales = as.vector(unlist(sales_data)),
Scenario = factor(rep(paste("Scenario", 1:num_scenarios), each = num_stores_per_scenario * num_months))
)
library(car)
leveneTest(Sales ~ Scenario, data = data)
by(data$Sales, data$Scenario, shapiro.test)
## data$Scenario: Scenario 1
##
## Shapiro-Wilk normality test
##
## data: dd[x, ]
## W = 0.95805, p-value = 0.1871
##
## ------------------------------------------------------------
## data$Scenario: Scenario 2
##
## Shapiro-Wilk normality test
##
## data: dd[x, ]
## W = 0.96123, p-value = 0.2345
##
## ------------------------------------------------------------
## data$Scenario: Scenario 3
##
## Shapiro-Wilk normality test
##
## data: dd[x, ]
## W = 0.97137, p-value = 0.464
anova_result <- aov(Sales ~ Scenario, data = data)
summary(anova_result)
## Df Sum Sq Mean Sq F value Pr(>F)
## Scenario 2 4340 2170.0 3.577 0.0314 *
## Residuals 105 63698 606.7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
| Df | Sum Sq | Mean Sq | F Value | Pr(>F) | |
|---|---|---|---|---|---|
| Scenario | 2 | 4340 | 2170.0 | 3.577 | 0.0314 |
| Residuals | 105 | 63698 | 606.7 |
When the sd was updated to 25 from 30, the p-value changed from 0.0554(5.54%) to 0.0314 (3.14%). The p-value being unique indicates that there is some differences in the sales, as the value has now dipped below the 5% common significance level. Furthermore this also suggests that promotional stratergies can perhaps have some impact on sales.
The report must be uploaded to Assignment 1 section in Canvas as a
PDF document with R codes and outputs showing. The
easiest way to do this is to:
1) Run all R chunks
2)
Preview your notebook in HTML (by
clicking Preview Notebook)
3) Open in Browser
(Chrome)
4) Right Click on the report in
Chrome
5) Click Print and Select the
Destination Option to Save as PDF.
6) Now upload
this PDF report as one single file via the Assignment 1 page in Canvas.
Remember to DELETE the instructional text provided in
the template. Failure to do this will INCREASE the SIMILARITY INDEX
reported in TURNITIN
If you have any questions regarding the assignment instructions or the R Markdown template, please post them on the discussion board.