Lab session 2: Solutions
Lab session 2: Solutions
Copyright By PowCoder代写 加微信 powcoder
RMIT University-MATH1302
23/03/2022
# List of packages required for this analysis
pkg <- c("mosaic", "ggplot2", "tidyverse", "car", "randtests", "stats", "dplyr", "agricolae")
# Check if packages are not installed and assign the
# names of the packages not installed to the variable new.pkg
new.pkg <- pkg[!(pkg %in% installed.packages())]
# If there are any packages in the list that aren't installed,
# install them
if (length(new.pkg)) {
install.packages(new.pkg, repos = "http://cran.rstudio.com")
# Load the packages into R
library(mosaic) # favstats()
library(ggplot2)
library(tidyverse) # summarise()
library(car) #ncvTest()
library(randtests) # runs.test()
library(stats) # TukeyHSD()
library(dplyr) # mutate()
library(agricolae) # LSD.test()
Question 1:
The effect of three different lubricating oils on fuel economy in diesel truck engines is being studied. Fuel economy is measured using brake-specific fuel consumption after the engine has been running for 15 minutes. Five different truck engines are available for the study, and the experimenters conduct the following randomized complete block design.
a) Create a data frame called “Fuel” in R to store the above data. Use str() command to display the structure of the data.
Insert the data frame.
Fuel<- data.frame(Oil = c(1,1,1,1,1,
2,2,2,2,2,
3,3,3,3,3),
Truck = c(1,2,3,4,5,
1,2,3,4,5,
1,2,3,4,5),
observations = c(0.500, 0.634, 0.487, 0.329, 0.512,
0.535, 0.675, 0.520, 0.435, 0.540,
0.513, 0.595, 0.488, 0.400, 0.510))
head(Fuel)
## Oil Truck observations
## 1 1 1 0.500
## 2 1 2 0.634
## 3 1 3 0.487
## 4 1 4 0.329
## 5 1 5 0.512
## 6 2 1 0.535
Export the data to the desktop
write.csv(Fuel,"/Users/abdulrahmanalamri/Desktop/Data/Fuel.csv", row.names = TRUE) # to Export DataFrame as a CSV file from RStudio to desktop
Display the structure of the data.
## 'data.frame': 15 obs. of 3 variables:
## $ Oil : num 1 1 1 1 1 2 2 2 2 2 ...
## $ Truck : num 1 2 3 4 5 1 2 3 4 5 ...
## $ observations: num 0.5 0.634 0.487 0.329 0.512 0.535 0.675 0.52 0.435 0.54 ...
b) Analyse the data from this experiment. State your hypotheses (use α = 5%) and draw conclusions.
First: we should treat the Oil and Teuck variables as factors.
Fuel <- Fuel %>% mutate(Oil = factor(Oil),
Truck = factor(Truck))
## ‘data.frame’: 15 obs. of 3 variables:
## $ Oil : Factor w/ 3 levels “1”,”2″,”3″: 1 1 1 1 1 2 2 2 2 2 …
## $ Truck : Factor w/ 5 levels “1”,”2″,”3″,”4″,..: 1 2 3 4 5 1 2 3 4 5 …
## $ observations: num 0.5 0.634 0.487 0.329 0.512 0.535 0.675 0.52 0.435 0.54 …
Second: state your hypothesis testing for the lubricating oils:
H0:τ1=τ2=τ3=0H0:τ1=τ2=τ3=0.
Ha:atleastoneτi≠0Ha:atleastoneτi≠0
where τiτi are the treatment effects of the three different lubricating oils.
Third: we might looks into the Descriptive Statistics:
Compute the mean and the SD (standard deviation) of the reaction time by different Ingredients:
group_by(Oil) %>%
summarise(n = n(),
mean = mean(observations, na.rm = TRUE),
sd = sd(observations, na.rm = TRUE),
stderr = sd/sqrt(n),
LCL = mean – qt(1 – (0.05 / 2), n – 1) * stderr,
UCL = mean + qt(1 – (0.05 / 2), n – 1) * stderr,
median = median(observations, na.rm = TRUE),
min = min(observations, na.rm = TRUE),
max = max(observations, na.rm = TRUE),
IQR = IQR(observations, na.rm = TRUE))
## # A tibble: 3 x 11
## Oil n mean sd stderr LCL UCL median min max IQR
##
## 1 1 5 0.492 0.109 0.0486 0.357 0.627 0.5 0.329 0.634 0.0250
## 2 2 5 0.541 0.0861 0.0385 0.434 0.648 0.535 0.435 0.675 0.0200
## 3 3 5 0.501 0.0697 0.0312 0.415 0.588 0.51 0.4 0.595 0.0250
Also, you can graphically check if there is a difference between the distribution using Box plot. Note that you include the jittered dot.
ggplot(Fuel, aes(x = Oil, y = observations, fill = Oil)) +
geom_boxplot() +
geom_jitter(shape = 15,
color = “steelblue”,
position = position_jitter(0.21)) +
theme_classic() +
labs( title= “Comparative Box Plot”)
Apply Complete Block design into the linear regression model.
model_1 <- lm(formula = observations ~ Oil + Truck, data=Fuel)
anova(model_1)
## Analysis of Variance Table
## Response: observations
## Df Sum Sq Mean Sq F value Pr(>F)
## Oil 2 0.006706 0.0033529 6.3527 0.02229 *
## Truck 4 0.092100 0.0230249 43.6257 1.781e-05 ***
## Residuals 8 0.004222 0.0005278
## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ‘ 1
Here, we might use two different approach in order to conclude on the output.
First: F-test.
H0H0 is rejected if F0>Fα,p−1,(p−2)(p−1)F0>Fα,p−1,(p−2)(p−1)
⸪F0=MSR/MSRes=6.3527⸪F0=MSR/MSRes=6.3527, We assumed α=0.05α=0.05
qf(p=.05, df1=2, df2=8, lower.tail=FALSE) # To find F critical value
## [1] 4.45897
⸪Fα,p−1,(p−2)(p−1)=F0.05,6,15=4.45897⸪Fα,p−1,(p−2)(p−1)=F0.05,6,15=4.45897
⸫F0≥Fα,p−1,(p−2)(p−1)⸫F0≥Fα,p−1,(p−2)(p−1)
As F0≥Fα,p−1,(p−2)(p−1)F0≥Fα,p−1,(p−2)(p−1), we reject the null hypothesis which means that at least one of three different lubricating oils effects is significant.
Second: p-value.
As the ANOVA table shown, the p-value of Ingredients is 0.022290.02229.So, we should reject the null hypothesis which mean the Oil factor is significant.
To obtain the estimates of the treatment effects we sould use Summary.
summary(model_1)
## lm(formula = observations ~ Oil + Truck, data = Fuel)
## Residuals:
## Min 1Q Median 3Q Max
## -0.039867 -0.008967 0.003133 0.010667 0.022333
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.496867 0.015694 31.660 1.08e-09 ***
## Oil2 0.048600 0.014530 3.345 0.010157 *
## Oil3 0.008800 0.014530 0.606 0.561529
## Truck2 0.118667 0.018758 6.326 0.000226 ***
## Truck3 -0.017667 0.018758 -0.942 0.373847
## Truck4 -0.128000 0.018758 -6.824 0.000135 ***
## Truck5 0.004667 0.018758 0.249 0.809795
## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ‘ 1
## Residual standard error: 0.02297 on 8 degrees of freedom
## Multiple R-squared: 0.959, Adjusted R-squared: 0.9283
## F-statistic: 31.2 on 6 and 8 DF, p-value: 3.958e-05
τ1=μ=0.496867,τ2=0.048600,τ3=0.008800τ1=μ=0.496867,τ2=0.048600,τ3=0.008800.
c) Construct the confidence intervals for each factor (for each level). Comment on the significance of the effect of each level based on this interval.
confint(model_1)
## 2.5 % 97.5 %
## (Intercept) 0.46067644 0.53305689
## Oil2 0.01509436 0.08210564
## Oil3 -0.02470564 0.04230564
## Truck2 0.07541107 0.16192226
## Truck3 -0.06092226 0.02558893
## Truck4 -0.17125559 -0.08474441
## Truck5 -0.03858893 0.04792226
We can check whether 00 is contained in the 95%95% confidence interval of all the estimators. Since the 95%95% confidence interval of Oil2 dose not contain 00, we do reject the null hypothesis at the significance level of 0.050.05, which mean that different lubricating oils will have a significant effect on fuel economy in diesel truck engines
d) Use the Fisher LSD method to make comparisons among the three lubricating oils to determine specifically which oils differ in brake-specific fuel consumption.
Use Fisher’s test to perform pairwise comparison.
LSD.test(y = model_1,
trt = “Oil”,
DFerror = model$df.residual,
MSerror = deviance(model)/model$df.residual,
alpha = 0.05,
group = TRUE,
console = TRUE)
## Study: model_1 ~ “Oil”
## LSD t Test for observations
## Mean Square Error: 0.0005277833
## Oil, means and individual ( 95 %) CI
## observations std r LCL UCL Min Max
## 1 0.4924 0.10865220 5 0.4687079 0.5160921 0.329 0.634
## 2 0.5410 0.08612491 5 0.5173079 0.5646921 0.435 0.675
## 3 0.5012 0.06969720 5 0.4775079 0.5248921 0.400 0.595
## Alpha: 0.05 ; DF Error: 8
## Critical Value of t: 2.306004
## least Significant Difference: 0.03350564
## Treatments with the same letter are not significantly different.
## observations groups
## 2 0.5410 a
## 3 0.5012 b
## 1 0.4924 b
As the test outlines, means that any pairs who do not share a letter are significantly different. Based on Fisher’s LSD method, All Oils have significantly different means. A more concise table that shows each pair-wise comparison is also provided below:
LSD.test(y = model_1,
trt = “Oil”,
DFerror = model$df.residual,
MSerror = deviance(model)/model$df.residual,
alpha = 0.05,
group = FALSE,
console = TRUE)
## Study: model_1 ~ “Oil”
## LSD t Test for observations
## Mean Square Error: 0.0005277833
## Oil, means and individual ( 95 %) CI
## observations std r LCL UCL Min Max
## 1 0.4924 0.10865220 5 0.4687079 0.5160921 0.329 0.634
## 2 0.5410 0.08612491 5 0.5173079 0.5646921 0.435 0.675
## 3 0.5012 0.06969720 5 0.4775079 0.5248921 0.400 0.595
## Alpha: 0.05 ; DF Error: 8
## Critical Value of t: 2.306004
## Comparison between treatments means
## difference pvalue signif. LCL UCL
## 1 – 2 -0.0486 0.0102 * -0.08210564 -0.01509436
## 1 – 3 -0.0088 0.5615 -0.04230564 0.02470564
## 2 – 3 0.0398 0.0255 * 0.00629436 0.07330564
From these intervals we can deduce that Oil 22 is significantly different than machines 11 and 33.
The Fisher’s LSD method shows us that the pairs 1−21−2 and 2−32−3 have means that are significantly different.
e) Analyse the residuals from this experiment and comment on model adequacy.
par(mfrow=c(2,2))
plot(model_1)
The normal probability plot shows the residuals follow an approximately straight line and that there are no outliers, verifying that the residuals are normally distributed. The residuals versus fits plots display points that have no recognisable pattern that fall randomly above and below zero, suggesting that the residuals are randomly distributed and verifying that they have constant variance and are homoscedastic.
Check assumptions statistically:
Residual Test:
a. Equality of variances – homogeneity
H0H0: variances are equal
H1H1: at least one variance is different
ncvTest(model_1)
## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 1.77572, Df = 1, p = 0.18268
The corresponding p-value (0.182680.18268) is more than α=0.05α=0.05. Therefore, we fail to reject the null hypothesis which means the model has a constant variance or “homogeneity”.
b. Normality of Residual
H0H0: Errors are normally distributed
H1H1: Errors are not normally distributed
y = rstudent(model_1)
shapiro.test(y)
## Shapiro-Wilk normality test
## data: y
## W = 0.82186, p-value = 0.007113
The p-Value 0.0071130.007113 is less than the significant level α=0.05α=0.05. So, we conclude that we do reject the null hypothesis test. Therefore, there is statistically significant evident that the stochastic component of this model is not normally distributed.
c. Auto Correlation Function:
H0H0: Errors are uncorrelated
H1H1: Errors are correlated
acf(model_1$residuals)
The ACF plot confirms that there is no any lags exceeded the horizontal dashed line. Therefore, we fail to reject the null hypothesis which means we have non-autocorrelation on residual model.
durbinWatsonTest(model_1)
## lag Autocorrelation D-W Statistic p-value
## 1 -0.3423932 2.682435 0.45
## Alternative hypothesis: rho != 0
The p-Value 0.4820.482 is much greater than the significant level α=0.05α=0.05. So, we conclude not to reject the null hypothesis test. Therefore, there is statistically significant evident that this model is non autocorrelation.
d. Randomness
H0H0: Errors are random
H1H1: Errors are not random
runs.test(rstudent(model_1))
## Runs Test
## data: rstudent(model_1)
## statistic = 2.2254, runs = 12, n1 = 7, n2 = 7, n = 14, p-value =
## 0.02605
## alternative hypothesis: nonrandomness
The corresponding p-value 0.026050.02605 is less than the significant level α=0.05α=0.05. Therefore, we reject null hypothesis. which means the independence of this model was nonrandomness.
Question 2:
The effect of five different ingredients (A, B, C, D, E) on the reaction time of a chemical process is being studied. Each batch of new material is only large enough to permit five runs to be made. Furthermore, each run requires approximately 1 and 1/2 hours, so only five runs can be made in one day. The experimenter decides to run the experiment as a Latin square so that day and batch effects may be systematically controlled. She obtains the data that follow.
Solution 2:
The aim of this question is to investigate the effect of five different ingredients (A, B, C, D, E) on the reaction time of a chemical process. Note that a Latin Square design of order t=5t=5 has been used in this experiment since we have two blocking factors and each treatment appears exactly one time in each row or column.
Note: A Latin square having the Latin letters sorted in the first row and column is called Normal or standard Latin square.
Insert the data frame.
Chemical<- data.frame(Batch = c(1,1,1,1,1,
2,2,2,2,2,
3,3,3,3,3,
4,4,4,4,4,
5,5,5,5,5), # rows
Day = c(1,2,3,4,5,
1,2,3,4,5,
1,2,3,4,5,
1,2,3,4,5,
1,2,3,4,5), # columns
Ingredients = c("A","B","D","C","E",
"C","E","A","D","B",
"B","A","C","E","D",
"D","C","E","B","A",
"E","D","B","A","C"), # Latin Letter
Time = c(8,7,1,7,3,
11,2,7,3,8,
4,9,10,1,5,
6,8,6,6,10,
4,2,3,8,8))
head(Chemical)
## Batch Day Ingredients Time
## 1 1 1 A 8
## 2 1 2 B 7
## 3 1 3 D 1
## 4 1 4 C 7
## 5 1 5 E 3
## 6 2 1 C 11
a) Use the ANOVA table to decide whether the five treatments are different. State your hypotheses and draw conclusions.
First: we should treat the Batch, Day and Ingredients variables as factors.
Chemical <- Chemical %>% mutate(Batch = factor(Batch),
Day = factor(Day),
Ingredients = factor(Ingredients))
str(Chemical)
## ‘data.frame’: 25 obs. of 4 variables:
## $ Batch : Factor w/ 5 levels “1”,”2″,”3″,”4″,..: 1 1 1 1 1 2 2 2 2 2 …
## $ Day : Factor w/ 5 levels “1”,”2″,”3″,”4″,..: 1 2 3 4 5 1 2 3 4 5 …
## $ Ingredients: Factor w/ 5 levels “A”,”B”,”C”,”D”,..: 1 2 4 3 5 3 5 1 4 2 …
## $ Time : num 8 7 1 7 3 11 2 7 3 8 …
Second: state your hypothesis testing for treatments.
H0:τ1=τ2=τ3=τ4=τ5=0H0:τ1=τ2=τ3=τ4=τ5=0.
Ha:atleastoneτi≠0Ha:atleastoneτi≠0
where τiτi are the treatment effects of the five different Ingredients.
Third: we might looks into the Descriptive Statistics:
Compute the mean and the SD (standard deviation) of the reaction time by different Ingredients:
Chemical %>%
group_by(Ingredients) %>%
summarise(n = n(),
mean = mean(Time, na.rm = TRUE),
sd = sd(Time, na.rm = TRUE),
stderr = sd/sqrt(n),
LCL = mean – qt(1 – (0.05 / 2), n – 1) * stderr,
UCL = mean + qt(1 – (0.05 / 2), n – 1) * stderr,
median = median(Time, na.rm = TRUE),
min = min(Time, na.rm = TRUE),
max = max(Time, na.rm = TRUE),
IQR = IQR(Time, na.rm = TRUE))
## # A tibble: 5 x 11
## Ingredients n mean sd stderr LCL UCL median min max IQR
##
## 1 A 5 8.4 1.14 0.510 6.98 9.82 8 7 10 1
## 2 B 5 5.6 2.07 0.927 3.03 8.17 6 3 8 3
## 3 C 5 8.8 1.64 0.735 6.76 10.8 8 7 11 2
## 4 D 5 3.4 2.07 0.927 0.825 5.97 3 1 6 3
## 5 E 5 3.2 1.92 0.860 0.812 5.59 3 1 6 2
Also, you can graphically check if there is a difference between the distribution using Box plot. Note that you include the jittered dot.
ggplot(Chemical, aes(x = Ingredients, y = Time, fill = Ingredients)) +
geom_boxplot() +
geom_jitter(shape = 15,
color = “steelblue”,
position = position_jitter(0.21)) +
theme_classic() +
labs( title= “Comparative Box Plot”)
Apply Latin Square design into the linear regression model.
model_2 <- lm(formula = Time ~ Batch + Day + Ingredients, data=Chemical)
anova(model_2)
## Analysis of Variance Table
## Response: Time
## Df Sum Sq Mean Sq F value Pr(>F)
## Batch 4 15.44 3.860 1.2345 0.3476182
## Day 4 12.24 3.060 0.9787 0.4550143
## Ingredients 4 141.44 35.360 11.3092 0.0004877 ***
## Residuals 12 37.52 3.127
## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ‘ 1
Here, we might use two different approach in order to conclude on the output.
First: F-test.
H0H0 is rejected if F0>Fα,p−1,(p−2)(p−1)F0>Fα,p−1,(p−2)(p−1)
⸪F0=MSR/MSRes=11.3092⸪F0=MSR/MSRes=11.3092, We assumed α=0.05α=0.05
qf(p=.05, df1=4, df2=12, lower.tail=FALSE) # To find F critical value
## [1] 3.259167
⸪Fα,p−1,(p−2)(p−1)=F0.05,3,15=3.259167⸪Fα,p−1,(p−2)(p−1)=F0.05,3,15=3.259167
⸫F0≥Fα,p−1,(p−2)(p−1)⸫F0≥Fα,p−1,(p−2)(p−1)
As F0≥Fα,p−1,(p−2)(p−1)F0≥Fα,p−1,(p−2)(p−1), we reject the null hypothesis which means that at least one of the five Ingredient effects is significant.
Second: p-value.
As the ANOVA table shown, the p-value of Ingredients is 0.00048770.0004877.So, we should reject the null hypothesis which mean the Ingredient factor is significant.
To obtain the estimates of the treatment effects we sould use Summary.
summary(model_2)
## lm(formula = Time ~ Batch + Day + Ingredients, data = Chemical)
## Residuals:
## Min 1Q Median 3Q Max
## -2.24 -1.24 -0.24 0.96 2.36
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.440 1.275 6.619 2.47e-05 ***
## Batch2 1.000 1.118 0.894 0.388803
## Batch3 0.600 1.118 0.537 0.601409
## Batch4 2.000 1.118 1.788 0.098971 .
## Batch5 -0.200 1.118 -0.179 0.861049
## Day2 -1.000 1.118 -0.894 0.388803
## Day3 -1.200 1.118 -1.073 0.304364
## Day4 -1.600 1.118 -1.431 0.178037
## Day5 0.200 1.
程序代写 CS代考 加微信: powcoder QQ: 1823890830 Email: powcoder@163.com