Lab 10
Homework
Question: Do light and fertilizer affect plant nutrient uptake?
IV: Light - 2 factors, fertilizer - 3 factors.
DV: Nutrient Uptake (continuous numerical)
Running a 2-factor ANOVA with 6 total levels.
1| Data read
light <- read.csv("Data/light_factorialanova.csv")
str(light) #turn characters to factors so tapply will work## 'data.frame': 180 obs. of 4 variables:
## $ Light : chr "Low" "Low" "Low" "Low" ...
## $ Fertilizer : chr "None" "None" "None" "None" ...
## $ Growth : num 5.5 4.86 5.65 6.52 4.77 ...
## $ Nutrient_Uptake: num 3.17 3.44 4.75 4.61 3.98 ...light$Light <- as.factor(light$Light)
light$Fertilizer <- as.factor(light$Fertilizer)
str(light) #great## 'data.frame': 180 obs. of 4 variables:
## $ Light : Factor w/ 2 levels "High","Low": 2 2 2 2 2 2 2 2 2 2 ...
## $ Fertilizer : Factor w/ 3 levels "None","Organic",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Growth : num 5.5 4.86 5.65 6.52 4.77 ...
## $ Nutrient_Uptake: num 3.17 3.44 4.75 4.61 3.98 ...2| Model fit
nutrientm <- lm(Nutrient_Uptake~Light*Fertilizer, data=light)
#check assumptions for residuals
par(mfrow=c(2,2))
plot(nutrientm)
#normality of residuals
shapiro.test(residuals(nutrientm)) #normality assumption met##
## Shapiro-Wilk normality test
##
## data: residuals(nutrientm)
## W = 0.99519, p-value = 0.8331#equal variance across species-sex groups
leveneTest(Nutrient_Uptake ~ interaction(Light, Fertilizer),data = light)#variance assumption met## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 5 2.1714 0.05941 .
## 174
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13| Summary & Anova
summary(nutrientm) #R2adj: 91% ##
## Call:
## lm(formula = Nutrient_Uptake ~ Light * Fertilizer, data = light)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.3521 -0.6284 -0.0091 0.5495 3.0323
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 14.0288 0.1777 78.967 < 2e-16 ***
## LightLow -9.7469 0.2512 -38.795 < 2e-16 ***
## FertilizerOrganic -2.0861 0.2512 -8.303 2.73e-14 ***
## FertilizerSynthetic -4.9822 0.2512 -19.830 < 2e-16 ***
## LightLow:FertilizerOrganic 5.5977 0.3553 15.755 < 2e-16 ***
## LightLow:FertilizerSynthetic 11.9285 0.3553 33.572 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.973 on 174 degrees of freedom
## Multiple R-squared: 0.9155, Adjusted R-squared: 0.9131
## F-statistic: 377.2 on 5 and 174 DF, p-value: < 2.2e-16Anova(nutrientm) #interaction term is significant## Anova Table (Type II tests)
##
## Response: Nutrient_Uptake
## Sum Sq Df F value Pr(>F)
## Light 686.14 1 724.683 < 2.2e-16 ***
## Fertilizer 30.90 2 16.317 3.201e-07 ***
## Light:Fertilizer 1068.51 2 564.264 < 2.2e-16 ***
## Residuals 164.75 174
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#NOTE if interaction term is not significant in this step, go back and take it out of model and just test main effects4| aov & TukeyHSD
nutrientaov <- aov(Nutrient_Uptake~Light*Fertilizer, data=light)
tuk <- TukeyHSD(nutrientaov)
tuk## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Nutrient_Uptake ~ Light * Fertilizer, data = light)
##
## $Light
## diff lwr upr p adj
## Low-High -3.904821 -4.191111 -3.618531 0
##
## $Fertilizer
## diff lwr upr p adj
## Organic-None 0.7127126 0.2927510 1.1326743 0.0002625
## Synthetic-None 0.9820468 0.5620851 1.4020084 0.0000003
## Synthetic-Organic 0.2693341 -0.1506275 0.6892958 0.2858761
##
## $`Light:Fertilizer`
## diff lwr upr p adj
## Low:None-High:None -9.7468825 -10.4708896 -9.022875495 0.0000000
## High:Organic-High:None -2.0861369 -2.8101440 -1.362129892 0.0000000
## Low:Organic-High:None -6.2353203 -6.9593274 -5.511313252 0.0000000
## High:Synthetic-High:None -4.9821965 -5.7062036 -4.258189477 0.0000000
## Low:Synthetic-High:None -2.8005925 -3.5245996 -2.076585443 0.0000000
## High:Organic-Low:None 7.6607456 6.9367385 8.384752658 0.0000000
## Low:Organic-Low:None 3.5115622 2.7875552 4.235569298 0.0000000
## High:Synthetic-Low:None 4.7646860 4.0406790 5.488693073 0.0000000
## Low:Synthetic-Low:None 6.9462901 6.2222830 7.670297107 0.0000000
## Low:Organic-High:Organic -4.1491834 -4.8731904 -3.425176305 0.0000000
## High:Synthetic-High:Organic -2.8960596 -3.6200666 -2.172052530 0.0000000
## Low:Synthetic-High:Organic -0.7144556 -1.4384626 0.009551504 0.0553886
## High:Synthetic-Low:Organic 1.2531238 0.5291167 1.977130830 0.0000216
## Low:Synthetic-Low:Organic 3.4347278 2.7107208 4.158734864 0.0000000
## Low:Synthetic-High:Synthetic 2.1816040 1.4575970 2.905611089 0.0000000letters <- multcompView::multcompLetters(tuk[["Light:Fertilizer"]][, "p adj"])
letters$Letters## Low:None High:Organic Low:Organic High:Synthetic Low:Synthetic
## "a" "b" "c" "d" "b"
## High:None
## "e"5| ggplot
#add significance letter
letters <- data.frame(Fertilizer = rep(c("None", "Organic", "Synthetic"), each = 2), Light = rep(c("High", "Low"), times = 3), label = c("a", "e", "b", "d", "c", "b"), height = c(16.5, 7, 14.7, 10, 13, 14.2))
#plot that is not box
ggplot(light, aes(Fertilizer,Nutrient_Uptake , fill=Light))+
geom_violin(position = position_dodge(width = 0.8))+
geom_text(data = letters, aes(Fertilizer, height, label = label), position = position_dodge(width = 0.8), show.legend = F) +
labs(x = "Fertilizer type", y = "Nutrient Uptake (mg/plant)", fill = "Light Level") +
theme_bw()+
theme(legend.position=c(0.92, 0.12), legend.background = element_rect(color="black"))+
scale_fill_manual(values = c("#488073ff", "#ceb526ff"))
Figure 1. Nutrient uptake (mg/plant) by fertilizer type and light level. The effect of fertilizer on nutrient uptake depended on light level, \(F_{2, 174}\) = 564.26, p < 0.01. Under high light, nutrient uptake was highest without fertilizer, whereas under low light, uptake was highest with synthetic fertilizer. The model explained 91.6% of the variation (\(R^2\) = .916). Different letters indicate significant differences among fertilizer–light groups (Tukey’s HSD, \(\alpha\) = .05).
Used to help create figure caption:
summary(nutrientm)$r.squared #multiply by 100 for percent explained## [1] 0.9155276Anova(nutrientm) #F, numerator df, and p; denominator df comes from Residuals## Anova Table (Type II tests)
##
## Response: Nutrient_Uptake
## Sum Sq Df F value Pr(>F)
## Light 686.14 1 724.683 < 2.2e-16 ***
## Fertilizer 30.90 2 16.317 3.201e-07 ***
## Light:Fertilizer 1068.51 2 564.264 < 2.2e-16 ***
## Residuals 164.75 174
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1This sentence is only here to practice pulling code from code chunck and putting it directly into my html paragraphs. The model explained 91.6% of the variation.
The first datapoint in the light dataframe has the following value for nutrient uptake: 3.171005 mg/plant.
Video Notes
Question: Do light and fertilizer affect plant growth? IV: Light - 2 factors, fertilizer - 3 factors. DV: Growth (continuous numerical) Running a 2-factor ANOVA with 6 total levels.
1| Data read
light <- read.csv("Data/light_factorialanova.csv")
str(light) #turn characters to factors so tapply will work## 'data.frame': 180 obs. of 4 variables:
## $ Light : chr "Low" "Low" "Low" "Low" ...
## $ Fertilizer : chr "None" "None" "None" "None" ...
## $ Growth : num 5.5 4.86 5.65 6.52 4.77 ...
## $ Nutrient_Uptake: num 3.17 3.44 4.75 4.61 3.98 ...light$Light <- as.factor(light$Light)
light$Fertilizer <- as.factor(light$Fertilizer)
str(light) #great## 'data.frame': 180 obs. of 4 variables:
## $ Light : Factor w/ 2 levels "High","Low": 2 2 2 2 2 2 2 2 2 2 ...
## $ Fertilizer : Factor w/ 3 levels "None","Organic",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Growth : num 5.5 4.86 5.65 6.52 4.77 ...
## $ Nutrient_Uptake: num 3.17 3.44 4.75 4.61 3.98 ...2| Model fit
growthm <- lm(Growth~Light*Fertilizer, data=light)
#check assumptions for residuals
par(mfrow=c(2,2))
plot(growthm)
#normality of residuals
shapiro.test(residuals(growthm)) #normality assumption met##
## Shapiro-Wilk normality test
##
## data: residuals(growthm)
## W = 0.99492, p-value = 0.8007#equal variance across species-sex groups
leveneTest(Growth ~ interaction(Light, Fertilizer),data = light)#variance assumption met## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 5 0.5499 0.7383
## 1743| Summary & Anova
summary(growthm) #R2adj: 96.1% ##
## Call:
## lm(formula = Growth ~ Light * Fertilizer, data = light)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.63263 -0.61351 -0.00418 0.55226 2.48349
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 14.9797 0.1737 86.220 < 2e-16 ***
## LightLow -10.1679 0.2457 -41.383 < 2e-16 ***
## FertilizerOrganic -3.0747 0.2457 -12.514 < 2e-16 ***
## FertilizerSynthetic 5.2990 0.2457 21.567 < 2e-16 ***
## LightLow:FertilizerOrganic 8.1417 0.3475 23.431 < 2e-16 ***
## LightLow:FertilizerSynthetic 1.9020 0.3475 5.474 1.52e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9516 on 174 degrees of freedom
## Multiple R-squared: 0.962, Adjusted R-squared: 0.9609
## F-statistic: 881.3 on 5 and 174 DF, p-value: < 2.2e-16Anova(growthm) #interaction term is significant## Anova Table (Type II tests)
##
## Response: Growth
## Sum Sq Df F value Pr(>F)
## Light 2093.06 1 2311.34 < 2.2e-16 ***
## Fertilizer 1353.17 2 747.15 < 2.2e-16 ***
## Light:Fertilizer 544.19 2 300.47 < 2.2e-16 ***
## Residuals 157.57 174
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#NOTE if interaction term is not significant in this step, go back and take it out of model and just test main effectsaov & TukeyHSD
growthaov <- aov(Growth~Light*Fertilizer, data=light)
TukeyHSD(growthaov)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Growth ~ Light * Fertilizer, data = light)
##
## $Light
## diff lwr upr p adj
## Low-High -6.82 -7.099983 -6.540018 0
##
## $Fertilizer
## diff lwr upr p adj
## Organic-None 0.996139 0.5854293 1.406849 1e-07
## Synthetic-None 6.250035 5.8393252 6.660745 0e+00
## Synthetic-Organic 5.253896 4.8431863 5.664606 0e+00
##
## $`Light:Fertilizer`
## diff lwr upr p adj
## Low:None-High:None -10.1678953 -10.875952 -9.4598384 0.000000
## High:Organic-High:None -3.0747065 -3.782763 -2.3666496 0.000000
## Low:Organic-High:None -5.1009108 -5.808968 -4.3928540 0.000000
## High:Synthetic-High:None 5.2990382 4.590981 6.0070951 0.000000
## Low:Synthetic-High:None -2.9668636 -3.674920 -2.2588067 0.000000
## High:Organic-Low:None 7.0931888 6.385132 7.8012457 0.000000
## Low:Organic-Low:None 5.0669844 4.358928 5.7750413 0.000000
## High:Synthetic-Low:None 15.4669335 14.758877 16.1749904 0.000000
## Low:Synthetic-Low:None 7.2010317 6.492975 7.9090886 0.000000
## Low:Organic-High:Organic -2.0262044 -2.734261 -1.3181475 0.000000
## High:Synthetic-High:Organic 8.3737447 7.665688 9.0818016 0.000000
## Low:Synthetic-High:Organic 0.1078429 -0.600214 0.8158998 0.997911
## High:Synthetic-Low:Organic 10.3999491 9.691892 11.1080060 0.000000
## Low:Synthetic-Low:Organic 2.1340472 1.425990 2.8421041 0.000000
## Low:Synthetic-High:Synthetic -8.2659018 -8.973959 -7.5578450 0.000000ggplot
#add significance letter
letters <- data.frame(Fertilizer = rep(c("None", "Organic", "Synthetic"), each = 2), Light = rep(c("High", "Low"), times = 3), label = c("b", "e", "c", "d", "a", "c"), height = c(18.5, 7.5, 15, 12.5, 23.5, 14.5))
#plot
ggplot(light, aes(Fertilizer,Growth , color=Light))+
geom_boxplot()+
geom_text(data = letters, aes(Fertilizer, height, label = label, color = Light), position = position_dodge(width = 0.8), show.legend = F) +
labs(x = "Fertilizer type", y = "Plant growth (cm)", color = "Light Level") +
theme_bw()+
theme(legend.position=c(0.92, 0.12), legend.background = element_rect(color="black"))
Model Table
| Analysis | Predictors | Example R formula | Where I used it |
|---|---|---|---|
| One-way / one-factor ANOVA | 1 categorical factor | Y ~ A |
Lab 09: deer diets, soil types, and plant disturbance treatments |
| Two-way / two-factor factorial ANOVA | 2 categorical factors | Y ~ A * B |
HW 10, Part A: millipede species × sex |
| Three-way / three-factor factorial ANOVA | 3 categorical factors | Y ~ A * B * C |
— |
| Four-way / four-factor factorial ANOVA | 4 categorical factors | Y ~ A * B * C * D |
— |
| n-way / n-factor factorial ANOVA | n categorical factors | Includes n factors | — |
| ANCOVA with a common slope | Categorical factor + numerical covariate | Y ~ A + X |
HW 10, Part B: mole-rat treatment + temperature |
| ANCOVA with an interaction | Categorical factor + numerical covariate + their interaction | Y ~ A * X |
HW 10, Part C: burn status × water depth predicting fish density |
| ANOVA with a block | Treatment factor + categorical block | Y ~ A + Block |
HW 10, Part B: mole-rat treatment + experimental run |
Factors are categorical predictors; levels are the categories within a factor. Numeric labels do not make a factor numerical.
In ANCOVA, the covariate is numerical.
A + X assumes a common slope across groups;
A * X allows slopes to differ.
In R formulas, A * B expands to
A + B + A:B, including both main effects and their
interaction.
| Question | Fixed block: + Block |
Random intercept: + (1 \| Block) |
|---|---|---|
| Example model | lm(SOC ~ Treatment + Block, data = soil) |
lmer(SOC ~ Treatment + (1 \| Block), data = soil) |
| Basic idea | Account for differences among these particular blocks. | Account for block differences using a shared distribution of block deviations. |
| Can blocks have different baselines? | Yes. | Yes. |
| What is a baseline? | Predicted SOC under the reference treatment, such as Control. | Predicted SOC under the reference treatment, such as Control. |
| What does the intercept mean? | Baseline for the reference block under Control. | Shared baseline under Control; individual blocks deviate from it. |
| What is a block estimate? | A difference from the reference block. | A deviation from the shared baseline. |
| Example | Intercept = 20 and BlockB = +6 predicts Control SOC of 26 in Block B. | Shared intercept = 20 and Block B deviation = +4 predicts Control SOC of 24 in Block B. |
| Where are individual block estimates? | Coefficient rows in summary(model), such as
BlockB. |
ranef(model)$Block. |
What does summary() show about random blocks? |
Not applicable. | Estimated variance and SD among block intercepts under “Random effects.” |
| What is partial pooling? | Block coefficients are not shrunk toward a shared baseline. | Less certain block estimates are pulled more toward the shared baseline. |
| Does the treatment effect vary among blocks? | No, not in this formula. | No, not in this formula. That would require a random treatment slope. |
| Does this change the observed data? | No. | No. |
| Important detail | Keep Block categorical. | 1 means intercept; \| Block means it
varies among blocks. |
Examples are illustrative, not results from my data. With additional numerical predictors, the intercept also refers to those predictors being zero.