Assignment 2_Heng Yan Ying
Heng Yan Ying
2026-09-22
Load packages
library(ggplot2)
library(dplyr)
library(tidyr)
library(ggpp) # offset scatter plot on box plot
library(car)
library(emmeans) # statsQuestion 1: Blocking
Load and review data
newts <- read.csv("C:/Users/charl/OneDrive - Nanyang Technological University/ES Mods/ES3307 Experimental Design/Assignment 2/Assignment 2 data/Newts.csv", header = TRUE)
# Check data structure and summary
str(newts)## 'data.frame': 120 obs. of 6 variables:
## $ SVL : num 108.9 94.7 96.9 110.1 103.7 ...
## $ Crest : num 10.12 5.62 2 15.04 3.88 ...
## $ Pond : int 4 7 7 10 6 8 7 6 7 2 ...
## $ Date : int 27 37 28 24 9 12 32 16 15 25 ...
## $ LSVL : num 2.04 1.98 1.99 2.04 2.02 ...
## $ LCREST: num 1.005 0.75 0.302 1.177 0.589 ...summary(newts)## SVL Crest Pond Date
## Min. : 87.60 Min. : 1.320 Min. : 1.00 Min. : 1.00
## 1st Qu.: 92.90 1st Qu.: 3.765 1st Qu.: 4.00 1st Qu.:11.50
## Median : 99.10 Median : 5.900 Median : 7.00 Median :22.00
## Mean : 99.64 Mean : 7.439 Mean : 6.08 Mean :21.99
## 3rd Qu.:105.75 3rd Qu.: 9.065 3rd Qu.: 8.00 3rd Qu.:32.00
## Max. :112.40 Max. :25.940 Max. :10.00 Max. :45.00
## NA's :33 NA's :33 NA's :33 NA's :33
## LSVL LCREST
## Min. :1.942 Min. :0.1200
## 1st Qu.:1.968 1st Qu.:0.5755
## Median :1.996 Median :0.7710
## Mean :1.997 Mean :0.7792
## 3rd Qu.:2.024 3rd Qu.:0.9575
## Max. :2.051 Max. :1.4140
## NA's :33 NA's :33Q1A
Show graphically how the size of the dorsal crest (LCREST) varies with the body size of the newts (LSVL) and POND. (Hint: Check your data structure to see if your variables are correctly specified i.e. continuous vs categorical.)
# Remove NA rows
newts <- na.omit(newts)
str(newts)## 'data.frame': 87 obs. of 6 variables:
## $ SVL : num 108.9 94.7 96.9 110.1 103.7 ...
## $ Crest : num 10.12 5.62 2 15.04 3.88 ...
## $ Pond : int 4 7 7 10 6 8 7 6 7 2 ...
## $ Date : int 27 37 28 24 9 12 32 16 15 25 ...
## $ LSVL : num 2.04 1.98 1.99 2.04 2.02 ...
## $ LCREST: num 1.005 0.75 0.302 1.177 0.589 ...
## - attr(*, "na.action")= 'omit' Named int [1:33] 88 89 90 91 92 93 94 95 96 97 ...
## ..- attr(*, "names")= chr [1:33] "88" "89" "90" "91" ...# Convert Pond to categorical variable
newts$Pond <- as.factor(newts$Pond)
str(newts)## 'data.frame': 87 obs. of 6 variables:
## $ SVL : num 108.9 94.7 96.9 110.1 103.7 ...
## $ Crest : num 10.12 5.62 2 15.04 3.88 ...
## $ Pond : Factor w/ 10 levels "1","2","3","4",..: 4 7 7 10 6 8 7 6 7 2 ...
## $ Date : int 27 37 28 24 9 12 32 16 15 25 ...
## $ LSVL : num 2.04 1.98 1.99 2.04 2.02 ...
## $ LCREST: num 1.005 0.75 0.302 1.177 0.589 ...
## - attr(*, "na.action")= 'omit' Named int [1:33] 88 89 90 91 92 93 94 95 96 97 ...
## ..- attr(*, "names")= chr [1:33] "88" "89" "90" "91" ...# Plot a scatter plot of LCREST against LSVL
newts %>%
ggplot(aes(x = LSVL, y = LCREST)) +
geom_point() +
labs(x = "Log of snout-vent length (mm)",
y = "Log of dorsal crest height (mm)") +
geom_smooth (method = "lm", color = "black") +
theme_minimal()
Fig. 1 Scatterplot of log of dorsal crest height (mm) against log of
snout-vent length (mm) with its linear model and standard error
bounds
# Plot a raincloud plot of dorsal crest height to compare across ponds
newts %>%
ggplot() +
geom_boxplot(aes(x = Pond, y = LCREST, fill = Pond),
width = 0.15,
outlier.shape = NA, # don't plot outliers, duplicated
alpha = 0.9) +
geom_jitter(aes(x = Pond, y = LCREST, colour = Pond), # scatterplot
position = position_jitternudge(width = 0.1,
height = 0,
x = -0.2), # offset pts
alpha = 0.5,
size = 1.2)+
labs(x = "Pond", y = "Log of dorsal crest height (mm)") +
theme_minimal()
Fig. 2 Raincloud plot of log of dorsal crest height (mm) against
ponds
Q1B
Analyse this dataset to investigate if the size of the dorsal crest (LCREST) reflects the body size of the newts (LSVL).
A linear model will be run to determine the relationship between size of dorsal crest and body size of newts. Diagnostic plots will be run on the linear model to determine if it is a suitable fit for the data.
newts_lm <- lm(LCREST ~ LSVL, newts)
anova(newts_lm)## Analysis of Variance Table
##
## Response: LCREST
## Df Sum Sq Mean Sq F value Pr(>F)
## LSVL 1 2.3894 2.38939 45.808 1.579e-09 ***
## Residuals 85 4.4337 0.05216
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(newts_lm)##
## Call:
## lm(formula = LCREST ~ LSVL, data = newts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.54974 -0.15546 0.00396 0.16713 0.52165
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -9.3806 1.5013 -6.248 1.59e-08 ***
## LSVL 5.0870 0.7516 6.768 1.58e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2284 on 85 degrees of freedom
## Multiple R-squared: 0.3502, Adjusted R-squared: 0.3425
## F-statistic: 45.81 on 1 and 85 DF, p-value: 1.579e-09There is a statistically significant relationship between LCREST and LSVL where p-value is less than 0.05 (F1,85=45.8, p-value = <0.001). As LSVL increases, LCREST increases (Fig. 1).
# Diagnostic plots of linear model's fit
par(mfrow = c(2, 2))
plot(newts_lm)
Fig 3. Diagnostic plots of linear model fitted to size of dorsal crest
against body size of newts
The Residuals vs Fitted line is roughly horizontal with points randomly scattered around 0, suggesting a linear relationship between LCREST and LSVL.
The Q-Q residuals graph shows the points generally follow the dotted line, suggesting a normal distribution of residuals.
The Scale-Location plot shows the data randomly scatter around the red line, suggesting homoscedasticity.
The Residuals vs Leverage graph shows there are no points lying outside of Cook’s distance (not drawn as out of graph bounds), which means there are no influential outliers in the data set.
Q1C
Run an analysis with POND as a blocking factor. Does POND seem to matter in this case? Should you keep it in the model?
newts %>%
ggplot(aes(x = LSVL, y = LCREST, colour = Pond)) + # group by Pond
geom_point() +
facet_wrap(~ Pond) # split plots by Pond
labs(title = "Pond",
x = "Log of snout-vent length (mm)",
y = "Log of dorsal crest height (mm)") +
theme_minimal()## NULLFig. 4 Scatter plot of LCREST against LSVL with Pond as a blocking factor
We will run the linear model again with and without Pond as a blocking factor.
# With Pond as blocking factor
newts_lm_pond <- lm(LCREST ~ Pond + LSVL, data = newts)
# Test the effect Pond + LSVL has on LCREST (Type 1 SS)
anova(newts_lm_pond)## Analysis of Variance Table
##
## Response: LCREST
## Df Sum Sq Mean Sq F value Pr(>F)
## Pond 9 0.3252 0.03613 0.6549 0.7466
## LSVL 1 2.3048 2.30483 41.7751 8.84e-09 ***
## Residuals 76 4.1931 0.05517
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Without Pond as blocking factor
newts_lm_pond2 <- lm(LCREST ~ LSVL + Pond, data = newts)
# Test the effect LSVL + Pond has on LCREST (Type 1 SS)
anova(newts_lm_pond2)## Analysis of Variance Table
##
## Response: LCREST
## Df Sum Sq Mean Sq F value Pr(>F)
## LSVL 1 2.3894 2.38939 43.3078 5.35e-09 ***
## Pond 9 0.2406 0.02674 0.4846 0.8807
## Residuals 76 4.1931 0.05517
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Test for Type 2 SS
Anova(newts_lm_pond, type = "2")## Anova Table (Type II tests)
##
## Response: LCREST
## Sum Sq Df F value Pr(>F)
## Pond 0.2406 9 0.4846 0.8807
## LSVL 2.3048 1 41.7751 8.84e-09 ***
## Residuals 4.1931 76
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Pond does not seem to have a statistically significant effect on LCREST when it was placed first in the model (F9,76 = 0.655, p = 0.747) or after LSVL (F9,76 = 0.485, p = 0.881).
Even after testing for Type 2 ANOVA which is order-independent unlike Type 1, Pond had no statistically significant effect on LCREST (F9,76 = 0.485, p = 0.881), while LSVL showed a significant effect (F1,76 = 41.8, p = <0.001), accounting for the variance.
Hence, as Pond does not contribute a significant effect to LCREST it is not necessary to use it as a blocking factor as each pond (Fig. 4) shows similar trends to the general trend (Fig. 1).
Q1D
To account for temporal variation, run an analysis with DATE as a covariate in the model. (Hint: think about where DATE should be placed in the model.)
newts %>%
ggplot(aes(x = LSVL, y = LCREST, colour = Date)) + # group by Date
geom_point() +
labs(title = "Date",
x = "Log of snout-vent length (mm)",
y = "Log of dorsal crest height (mm)") +
geom_smooth (method = "lm", color = "black") +
theme_minimal()
Fig. 5 LCREST against LSVL organised by Date
# Place Date first in model
newts_lm_date = lm(LCREST ~ Date + Pond + LSVL, data = newts)
anova(newts_lm_date)## Analysis of Variance Table
##
## Response: LCREST
## Df Sum Sq Mean Sq F value Pr(>F)
## Date 1 0.0310 0.03101 0.5613 0.4561
## Pond 9 0.3223 0.03581 0.6481 0.7524
## LSVL 1 2.3257 2.32567 42.0896 8.336e-09 ***
## Residuals 75 4.1441 0.05526
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Place Pond first, LSVL last (since it explains most var)
newts_lm_date2 = lm(LCREST ~ Pond + Date + LSVL, data = newts)
anova(newts_lm_date2)## Analysis of Variance Table
##
## Response: LCREST
## Df Sum Sq Mean Sq F value Pr(>F)
## Pond 9 0.3252 0.03613 0.6539 0.7474
## Date 1 0.0281 0.02810 0.5086 0.4779
## LSVL 1 2.3257 2.32567 42.0896 8.336e-09 ***
## Residuals 75 4.1441 0.05526
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Test for Type 2 SS
Anova(newts_lm_date, type = "2")## Anova Table (Type II tests)
##
## Response: LCREST
## Sum Sq Df F value Pr(>F)
## Date 0.0490 1 0.8859 0.3496
## Pond 0.2222 9 0.4467 0.9050
## LSVL 2.3257 1 42.0896 8.336e-09 ***
## Residuals 4.1441 75
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Date does not seem to have a statistically significant effect on LCREST when it was placed first in the model (F1,75 = 0.561, p = 0.456) or second, after Pond (F1,75 = 0.509, p = 0.478).
Even after testing for Type 2 ANOVA which is order-independent unlike Type 1, Date had no statistically significant effect on LCREST (F1,75 = 0.886, p = 0.350). Thus, Date should not be included in the model as a covariate.
Question 2: Psuedoreplication
Load and review data
sheeps <- read.csv("C:/Users/charl/OneDrive - Nanyang Technological University/ES Mods/ES3307 Experimental Design/Assignment 2/Assignment 2 data/sheep.csv", header = TRUE)
# Check data structure and summary
str(sheeps)## 'data.frame': 120 obs. of 10 variables:
## $ Duration : int 39 40 39 34 31 37 38 36 36 37 ...
## $ nlookups : int 3 7 3 7 6 5 3 4 2 5 ...
## $ sex : int 1 1 1 1 1 1 1 1 1 1 ...
## $ sheep : int 1 1 1 1 1 1 1 1 1 1 ...
## $ sheepn : int 1 1 1 1 1 1 1 1 1 1 ...
## $ obsper : int 1 2 3 4 5 6 7 8 9 10 ...
## $ luprate : num 0.0769 0.175 0.0769 0.2059 0.1935 ...
## $ Sex2 : int 1 1 1 2 2 2 NA NA NA NA ...
## $ sheep2 : int 1 2 3 4 5 6 NA NA NA NA ...
## $ avluprate: num 0.133 0.235 0.168 0.195 0.217 ...summary(sheeps)## Duration nlookups sex sheep sheepn
## Min. :27.00 Min. : 2.00 Min. :1.0 Min. :1.0 Min. :1
## 1st Qu.:34.75 1st Qu.: 6.00 1st Qu.:1.0 1st Qu.:2.0 1st Qu.:1
## Median :38.00 Median : 7.00 Median :1.5 Median :3.5 Median :2
## Mean :38.09 Mean : 7.85 Mean :1.5 Mean :3.5 Mean :2
## 3rd Qu.:41.00 3rd Qu.:10.00 3rd Qu.:2.0 3rd Qu.:5.0 3rd Qu.:3
## Max. :53.00 Max. :14.00 Max. :2.0 Max. :6.0 Max. :3
##
## obsper luprate Sex2 sheep2
## Min. : 1.00 Min. :0.05556 Min. :1.0 Min. :1.00
## 1st Qu.: 5.75 1st Qu.:0.15790 1st Qu.:1.0 1st Qu.:2.25
## Median :10.50 Median :0.20521 Median :1.5 Median :3.50
## Mean :10.50 Mean :0.21164 Mean :1.5 Mean :3.50
## 3rd Qu.:15.25 3rd Qu.:0.25659 3rd Qu.:2.0 3rd Qu.:4.75
## Max. :20.00 Max. :0.43333 Max. :2.0 Max. :6.00
## NA's :114 NA's :114
## avluprate
## Min. :0.1327
## 1st Qu.:0.1744
## Median :0.2059
## Mean :0.2116
## 3rd Qu.:0.2304
## Max. :0.3228
## NA's :114# Clean data
sheeps <- sheeps[-c(5,8:10)] # Remove unwanted columns
sheeps <- sheeps %>% # Convert Sex into male and female
mutate(sex = case_when(sex == 1 ~ "Female",
sex == 2 ~ "Male"))
# Convert obsper and sheep to factor variables
sheeps$obsper <- as.factor(sheeps$obsper)
sheeps$sheep <- as.factor(sheeps$sheep)
str(sheeps)## 'data.frame': 120 obs. of 6 variables:
## $ Duration: int 39 40 39 34 31 37 38 36 36 37 ...
## $ nlookups: int 3 7 3 7 6 5 3 4 2 5 ...
## $ sex : chr "Female" "Female" "Female" "Female" ...
## $ sheep : Factor w/ 6 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ obsper : Factor w/ 20 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ luprate : num 0.0769 0.175 0.0769 0.2059 0.1935 ...Q2A
Fit a linear model of LUPRATE = OBSPER + SEX. Criticise this analysis.
# Split number of observations in half to plot 2 halves
sheeps <- sheeps %>%
mutate(obsper_half = ifelse(as.numeric(obsper) <= 10, "1-10", "11-20"))
sheeps %>%
ggplot() +
geom_boxplot(aes(x = obsper, y = luprate, fill = sex),
width = 0.35,
outlier.shape = NA, # don't plot outliers, duplicated
alpha = 0.9) +
# Split scatter points by gender
geom_jitter(data = filter(sheeps, sex == "Female"),
aes(x = obsper, y = luprate, colour = sex),
position = position_jitternudge(width = 0,
height = 0,
x = -0.25),
alpha = 0.5,
size = 0.9) +
geom_jitter(data = filter(sheeps, sex == "Male"),
aes(x = obsper, y = luprate, colour = sex),
position = position_jitternudge(width = 0,
height = 0,
x = 0.25),
alpha = 0.5,
size = 0.9) +
labs(x = "Observation Period", y = "Look-up rate") +
# Join both observation halves together
facet_wrap(~ obsper_half, ncol = 1, scales = "free_x") +
theme_minimal()
Fig. 6 Look-up rate against Observation Period for Female and Male
sheep
We will fit a linear model of LUPRATE = OBSPER + SEX.
sheeps_lm <- lm(luprate ~ obsper + sex, data = sheeps)
anova(sheeps_lm)## Analysis of Variance Table
##
## Response: luprate
## Df Sum Sq Mean Sq F value Pr(>F)
## obsper 19 0.19192 0.010101 2.4082 0.002655 **
## sex 1 0.13282 0.132816 31.6652 1.707e-07 ***
## Residuals 99 0.41524 0.004194
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1According to ANOVA, both observation period and sex had a significant effect on look-up rate as both p-value was less than 0.05 (observation period: F19,99 = 2.41, p = 0.00266; sex: F1,99 = 31.7, p = <0.001).
However, there is pseudoreplication in this experiment as the look-up rate was taken from the same 6 sheeps. This is shown in the high residual df of 99 even though there were only 6 study subjects.
Q2B
Suggest and execute a more appropriate analysis.
Instead of taking each observation as an independent data point which would result in pseudoreplication, each sheep should be taken as one independent observation. This can be done by averaging look-up rate for each sheep. Then, they would be grouped by sex to compare average look-up rate between male and female sheeps.
# Group by sheep and sex
sheeps_avg <- sheeps %>%
group_by(sheep, sex) %>%
summarize(luprate_avg = mean(luprate))
# Box plot comparing avg luprate btwn sex
sheeps_avg %>%
ggplot(aes(x = sex, y = luprate_avg)) +
geom_boxplot(width = 0.2, fill = c("pink","lightblue")) +
labs(x = "Sex", y = "Average Look-up Rate")
Fig. 7 Box plot of average look-up rate for female and male sheeps
As n = 3 is a small sample size, it is difficult to determine if the distribution is normal or not. Hence, the non-parametric Mann-Whitney Test will be used to compare the average look-up time between sexes.
wilcox.test(luprate_avg ~ sex, data = sheeps_avg)##
## Wilcoxon rank sum exact test
##
## data: luprate_avg by sex
## W = 2, p-value = 0.4
## alternative hypothesis: true location shift is not equal to 0From the Mann-Whitney test, the average look-up rate between sexes is not statistically significant (W = 2, p = 0.4) as the p-value is greater than 0.05.
Question 3: Sums of Squares
Load and review data
weight <- read.csv("C:/Users/charl/OneDrive - Nanyang Technological University/ES Mods/ES3307 Experimental Design/Assignment 2/Assignment 2 data/weight.csv", header = TRUE)
# Check data structure and summary
str(weight)## 'data.frame': 39 obs. of 3 variables:
## $ FOREARM: num 4.87 4.67 3 3.54 1.72 ...
## $ HT : num 159 164 158 148 155 ...
## $ WT : num 70.7 63.9 60.4 56.7 56.6 ...summary(weight)## FOREARM HT WT
## Min. :0.232 Min. :143.7 Min. :40.43
## 1st Qu.:3.477 1st Qu.:154.2 1st Qu.:59.37
## Median :4.868 Median :157.1 Median :63.55
## Mean :5.071 Mean :157.9 Mean :63.19
## 3rd Qu.:6.577 3rd Qu.:160.9 3rd Qu.:66.58
## Max. :9.811 Max. :171.2 Max. :80.93Q3A
Using two separate plots, show graphically how FOREARM varies with HT and WT.
weight %>%
ggplot(aes(x = HT, y = FOREARM)) +
geom_point() +
labs(x = "Height", y = "Forearm thickness") +
geom_smooth (method = "lm", color = "black") +
theme_minimal()
Fig. 8 Scatter plot of forearm thickness against height
weight %>%
ggplot(aes(x = WT, y = FOREARM)) +
geom_point() +
labs(x = "Weight", y = "Forearm thickness") +
geom_smooth (method = "lm", color = "black") +
theme_minimal()
Fig. 9 Scatter plot of forearm thickness against weight
Q3B
Taking FOREARM as the response variable, which of the two explanatory variables, HT and WT, is the better predictor of obesity when used alone in a linear model?
# ANOVA of HT as explanatory variable
weight_lm_HT <- lm(FOREARM ~ HT, data = weight)
anova(weight_lm_HT)## Analysis of Variance Table
##
## Response: FOREARM
## Df Sum Sq Mean Sq F value Pr(>F)
## HT 1 0.944 0.9439 0.1754 0.6778
## Residuals 37 199.094 5.3809# ANOVA of WT as explanatory variable
weight_lm_WT <- lm(FOREARM ~ WT, data = weight)
anova(weight_lm_WT)## Analysis of Variance Table
##
## Response: FOREARM
## Df Sum Sq Mean Sq F value Pr(>F)
## WT 1 59.137 59.137 15.529 0.000347 ***
## Residuals 37 140.901 3.808
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1When used in a linear model, weight is a better predictor of obesity as p-value is less than 0.05 (F1,37 = 15.5, p-value = 0.000347) while p-value for height as an explanatory variable is greater than 0.05 (F1,37 = 0.175, p-value = 0.678). Additionally, from Fig 9, there is a clear upward trend as weight increases forearm thickness increases, but from Fig 8 the trend between height and forearm thickness is not as clear.
Q3C
If both HT and WT are included as main effects (no interactions) in a linear model, do they increase or detract from each others’ informativeness, and why? (hint: try running with different sums of squares).
# Run a model with HT then WT
obesity_lm1 <- lm(FOREARM ~ HT + WT, data = weight)
anova(obesity_lm1)## Analysis of Variance Table
##
## Response: FOREARM
## Df Sum Sq Mean Sq F value Pr(>F)
## HT 1 0.944 0.944 0.2926 0.5919
## WT 1 82.970 82.970 25.7220 1.207e-05 ***
## Residuals 36 116.124 3.226
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Run a model with WT then HT
obesity_lm2 <- lm(FOREARM ~ WT + HT, data = weight)
anova(obesity_lm2)## Analysis of Variance Table
##
## Response: FOREARM
## Df Sum Sq Mean Sq F value Pr(>F)
## WT 1 59.137 59.137 18.3334 0.0001314 ***
## HT 1 24.777 24.777 7.6813 0.0087755 **
## Residuals 36 116.124 3.226
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Test for Type 2 ANOVA
Anova(obesity_lm1, type = "2")## Anova Table (Type II tests)
##
## Response: FOREARM
## Sum Sq Df F value Pr(>F)
## HT 24.777 1 7.6813 0.008775 **
## WT 82.970 1 25.7220 1.207e-05 ***
## Residuals 116.124 36
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1When putting HT first then WT, HT does not have a significant effect (F1,36 = 0.293, p = 0.591) on obesity as p-value is greater than 0.05, while WT does (F1,36 = 25.7, p = <0.001) as p-value is less than 0.05.
However, when putting WT first then HT, both WT and HT turn out to have a significant effect on obesity (WT: F1,36 = 18.3, p = 0.000131; HT: F1,36 = 7.68, p = 0.00878) as p-value is less than 0.05.
When Type 2 ANOVA was run which is order-independent, both HT and WT have significant effect on obesity (HT: F1,36 = 7.68, p = 0.00878; WT: F1,36 = 25.7, p = <0.001).
Q3D
What happens when you put HT vs WT first in the models?
The order of the variable put into the model affects whether Type 1 ANOVA returns a significant effect on the response variable or not. When HT was put first, HT was not a significant effect but WT was. When WT was put first, both HT and WT were a significant effect on obesity.
Question 4: Orthogonality
Load and review data
blooms <- read.csv("C:/Users/charl/OneDrive - Nanyang Technological University/ES Mods/ES3307 Experimental Design/Assignment 2/Assignment 2 data/blooms.csv", header = TRUE)
# Check data structure and summary
str(blooms)## 'data.frame': 36 obs. of 9 variables:
## $ SQBLOOMS: num 4.36 3.32 3.61 4.12 4.47 ...
## $ BED : int 1 1 1 1 1 1 1 1 1 1 ...
## $ WATER : int 1 1 1 1 2 2 2 2 3 3 ...
## $ SHADE : int 1 2 3 4 1 2 3 4 1 2 ...
## $ X : logi NA NA NA NA NA NA ...
## $ SQ2 : num 4.36 3.32 3.61 4.12 4.47 ...
## $ B2 : int 1 1 1 1 1 1 1 1 1 1 ...
## $ W2 : int 1 1 1 1 2 2 2 2 3 3 ...
## $ S2 : int 1 2 3 4 1 2 3 4 1 2 ...summary(blooms)## SQBLOOMS BED WATER SHADE X
## Min. :2.449 Min. :1 Min. :1 Min. :1.00 Mode:logical
## 1st Qu.:3.571 1st Qu.:1 1st Qu.:1 1st Qu.:1.75 NA's:36
## Median :4.061 Median :2 Median :2 Median :2.50
## Mean :4.029 Mean :2 Mean :2 Mean :2.50
## 3rd Qu.:4.583 3rd Qu.:3 3rd Qu.:3 3rd Qu.:3.25
## Max. :5.385 Max. :3 Max. :3 Max. :4.00
##
## SQ2 B2 W2 S2
## Min. :2.828 Min. :1.00 Min. :1 Min. :1.000
## 1st Qu.:3.606 1st Qu.:1.00 1st Qu.:1 1st Qu.:1.000
## Median :4.123 Median :2.00 Median :2 Median :2.000
## Mean :4.094 Mean :1.97 Mean :2 Mean :2.424
## 3rd Qu.:4.583 3rd Qu.:3.00 3rd Qu.:3 3rd Qu.:3.000
## Max. :5.385 Max. :3.00 Max. :3 Max. :4.000
## NA's :3 NA's :3 NA's :3 NA's :3# Remove column X
blooms <- blooms[-c(5)]
# Convert BED, WATER, SHADE, B2, W2, S2 into factor variables
blooms$BED <- as.factor(blooms$BED)
blooms$WATER <- as.factor(blooms$WATER)
blooms$SHADE <- as.factor(blooms$SHADE)
blooms$B2 <- as.factor(blooms$B2)
blooms$W2 <- as.factor(blooms$W2)
blooms$S2 <- as.factor(blooms$S2)
# Split data into 2 tables
blooms_1 <- blooms[c(1:4)]
blooms_2 <- blooms[c(5:8)]
# Remove NA rows
blooms_2 <- na.omit(blooms_2)
str(blooms_1)## 'data.frame': 36 obs. of 4 variables:
## $ SQBLOOMS: num 4.36 3.32 3.61 4.12 4.47 ...
## $ BED : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
## $ WATER : Factor w/ 3 levels "1","2","3": 1 1 1 1 2 2 2 2 3 3 ...
## $ SHADE : Factor w/ 4 levels "1","2","3","4": 1 2 3 4 1 2 3 4 1 2 ...str(blooms_2)## 'data.frame': 33 obs. of 4 variables:
## $ SQ2: num 4.36 3.32 3.61 4.12 4.47 ...
## $ B2 : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
## $ W2 : Factor w/ 3 levels "1","2","3": 1 1 1 1 2 2 2 2 3 3 ...
## $ S2 : Factor w/ 4 levels "1","2","3","4": 1 2 3 4 1 2 3 4 1 2 ...
## - attr(*, "na.action")= 'omit' Named int [1:3] 34 35 36
## ..- attr(*, "names")= chr [1:3] "34" "35" "36"Q4A
Show graphically how the number of blooms varies with level of water and level of shade.
# Box plot of how blooms vary with level of water and shade
blooms %>%
ggplot() +
geom_boxplot(aes(x = WATER, y = SQBLOOMS, fill = SHADE),
width = 0.5,
outlier.shape = NA, # don't plot outliers, duplicated
alpha = 0.9) +
labs(x = "Level of Watering",
y = "Square root of blooms",
fill = "Shade Level") +
theme_minimal()
Fig. 10 Box plot of square root of blooms against level of watering and
separated by shade levels
Q4B
Fit a linear model of SQBLOOMS = BED + WATER + SHADE. Is the analysis orthogonal?
blooms_lm <- lm(SQBLOOMS ~ BED + WATER + SHADE, data = blooms_1)
anova(blooms_lm)## Analysis of Variance Table
##
## Response: SQBLOOMS
## Df Sum Sq Mean Sq F value Pr(>F)
## BED 2 4.1323 2.06614 9.4570 0.0007277 ***
## WATER 2 3.7153 1.85767 8.5029 0.0013016 **
## SHADE 3 1.6465 0.54882 2.5120 0.0789451 .
## Residuals 28 6.1173 0.21848
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Anova(blooms_lm, type = "2")## Anova Table (Type II tests)
##
## Response: SQBLOOMS
## Sum Sq Df F value Pr(>F)
## BED 4.1323 2 9.4570 0.0007277 ***
## WATER 3.7153 2 8.5029 0.0013016 **
## SHADE 1.6465 3 2.5120 0.0789451 .
## Residuals 6.1173 28
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Since Type 1 and Type 2 ANOVA returns the same F and p-value, the design is orthogonal and balanced.
Q4C
Fit a linear model of SQBLOOMS = WATER + SHADE. How were the residual error and sums of squares affected in this analysis? Was it worthwhile blocking for BED? If so, why?
blooms_lm2 <- lm(SQBLOOMS ~ WATER + SHADE, data = blooms_1)
anova(blooms_lm2)## Analysis of Variance Table
##
## Response: SQBLOOMS
## Df Sum Sq Mean Sq F value Pr(>F)
## WATER 2 3.7153 1.85767 5.4373 0.009661 **
## SHADE 3 1.6465 0.54882 1.6064 0.208609
## Residuals 30 10.2496 0.34165
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Anova(blooms_lm2)## Anova Table (Type II tests)
##
## Response: SQBLOOMS
## Sum Sq Df F value Pr(>F)
## WATER 3.7153 2 5.4373 0.009661 **
## SHADE 1.6465 3 1.6064 0.208609
## Residuals 10.2496 30
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The residuals sum of squares (SS) increased by 4.13 and residuals degrees of freedom (df) increased by 2. This corresponds exactly to the df and SS contributed by BED in the SQBLOOMS = BED + WATER + SHADE model. When BED was removed as a covariate, its df and SS got redistributed to the unexplained error. Hence, it was worthwhile blocking for BED as it is able to explain some of the unexplained error, which reduces variance from the residual error term.
Q4D
Fit the model SQ2 = B2 + W2 + S2. Which parts of the output are different in this third model, but were the same in the first two, and why?
blooms_2_lm <- lm(SQ2 ~ B2 + W2 + S2, data = blooms_2)
anova(blooms_2_lm)## Analysis of Variance Table
##
## Response: SQ2
## Df Sum Sq Mean Sq F value Pr(>F)
## B2 2 2.7626 1.38129 8.3790 0.001641 **
## W2 2 5.0793 2.53966 15.4057 4.367e-05 ***
## S2 3 0.8072 0.26906 1.6321 0.207156
## Residuals 25 4.1213 0.16485
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Anova(blooms_2_lm, type = "2")## Anova Table (Type II tests)
##
## Response: SQ2
## Sum Sq Df F value Pr(>F)
## B2 2.6490 2 8.0346 0.00202 **
## W2 4.6764 2 14.1836 7.644e-05 ***
## S2 0.8072 3 1.6321 0.20716
## Residuals 4.1213 25
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1For this model, when Type 1 and 2 ANOVA were run it did not return the same F and p values for the corresponding terms B2 and S2 which is attributed by less data points (33 observations) compared to the first 2 models which were orthogonal (36 observations). Since this design is not orthogonal, running both ANOVA types causes unbalanced comparison across BED and WATER combinations. Thus, this design will be affected by the ordering of which terms appear first in the linear model unlike the first 2.