growth.moo <- read.csv("/Users/aaronsi/Documents/Documents - Aaron’s MacBook Pro/4259A/data/growth.csv")
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(ggplot2)
library(ggfortify)
glimpse(growth.moo)
## Rows: 48
## Columns: 3
## $ supplement <chr> "supergain", "supergain", "supergain", "supergain", "contro…
## $ diet <chr> "wheat", "wheat", "wheat", "wheat", "wheat", "wheat", "whea…
## $ gain <dbl> 17.37125, 16.81489, 18.08184, 15.78175, 17.70656, 18.22717,…
growth.moo <- mutate(growth.moo,
diet = factor(diet),
supplement = factor(supplement))
levels(growth.moo$diet)
## [1] "barley" "oats" "wheat"
levels(growth.moo$supplement)
## [1] "agrimore" "control" "supergain" "supersupp"
growth.moo <-
mutate(growth.moo,
supplement = relevel(supplement, ref = "control"))
levels(growth.moo$supplement)
## [1] "control" "agrimore" "supergain" "supersupp"
sumMoo <- growth.moo %>%
group_by(diet, supplement) %>%
summarise(meanGrow = mean(gain))
## `summarise()` has regrouped the output.
## ℹ Summaries were computed grouped by diet and supplement.
## ℹ Output is grouped by diet.
## ℹ Use `summarise(.groups = "drop_last")` to silence this message.
## ℹ Use `summarise(.by = c(diet, supplement))` for per-operation grouping
## (`?dplyr::dplyr_by`) instead.
sumMoo
## # A tibble: 12 × 3
## # Groups: diet [3]
## diet supplement meanGrow
## <fct> <fct> <dbl>
## 1 barley control 23.3
## 2 barley agrimore 26.3
## 3 barley supergain 22.5
## 4 barley supersupp 25.6
## 5 oats control 20.5
## 6 oats agrimore 23.3
## 7 oats supergain 19.7
## 8 oats supersupp 21.9
## 9 wheat control 17.4
## 10 wheat agrimore 19.6
## 11 wheat supergain 17.0
## 12 wheat supersupp 19.7
ggplot(sumMoo, aes(x = supplement, y = meanGrow)) +
geom_point() +
theme_bw()
ggplot(sumMoo, aes(x = supplement, y = meanGrow,
colour = diet, group = diet)) +
geom_point() +
geom_line() +
theme_bw()
model_cow <- lm(gain ~ diet * supplement, data = growth.moo)
autoplot(model_cow, smooth.colour = NA)
## Warning: `fortify(<lm>)` was deprecated in ggplot2 4.0.0.
## ℹ Please use `broom::augment(<lm>)` instead.
## ℹ The deprecated feature was likely used in the ggfortify package.
## Please report the issue at <https://github.com/sinhrks/ggfortify/issues>.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## ℹ The deprecated feature was likely used in the ggfortify package.
## Please report the issue at <https://github.com/sinhrks/ggfortify/issues>.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: Removed 48 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 48 rows containing missing values or values outside the scale range
## (`geom_line()`).
anova(model_cow)
## Analysis of Variance Table
##
## Response: gain
## Df Sum Sq Mean Sq F value Pr(>F)
## diet 2 287.171 143.586 83.5201 2.999e-14 ***
## supplement 3 91.881 30.627 17.8150 2.952e-07 ***
## diet:supplement 6 3.406 0.568 0.3302 0.9166
## Residuals 36 61.890 1.719
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model_cow)
##
## Call:
## lm(formula = gain ~ diet * supplement, data = growth.moo)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.48756 -1.00368 -0.07452 1.03496 2.68069
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 23.2966499 0.6555863 35.536 < 2e-16 ***
## dietoats -2.8029851 0.9271390 -3.023 0.00459 **
## dietwheat -5.8911317 0.9271390 -6.354 2.34e-07 ***
## supplementagrimore 3.0518277 0.9271390 3.292 0.00224 **
## supplementsupergain -0.8305263 0.9271390 -0.896 0.37631
## supplementsupersupp 2.2786527 0.9271390 2.458 0.01893 *
## dietoats:supplementagrimore -0.2471088 1.3111726 -0.188 0.85157
## dietwheat:supplementagrimore -0.8182729 1.3111726 -0.624 0.53651
## dietoats:supplementsupergain -0.0001351 1.3111726 0.000 0.99992
## dietwheat:supplementsupergain 0.4374395 1.3111726 0.334 0.74060
## dietoats:supplementsupersupp -0.9120830 1.3111726 -0.696 0.49113
## dietwheat:supplementsupersupp -0.0158299 1.3111726 -0.012 0.99043
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.311 on 36 degrees of freedom
## Multiple R-squared: 0.8607, Adjusted R-squared: 0.8182
## F-statistic: 20.22 on 11 and 36 DF, p-value: 3.295e-12
sumMoo <- growth.moo %>%
group_by(diet, supplement) %>%
summarise(
meanGrow = mean(gain),
seGrow = sd(gain)/sqrt(n())
)
## `summarise()` has regrouped the output.
## ℹ Summaries were computed grouped by diet and supplement.
## ℹ Output is grouped by diet.
## ℹ Use `summarise(.groups = "drop_last")` to silence this message.
## ℹ Use `summarise(.by = c(diet, supplement))` for per-operation grouping
## (`?dplyr::dplyr_by`) instead.
ggplot(sumMoo, aes(x = supplement, y = meanGrow,
colour = diet, group = diet)) +
geom_point() +
geom_line() +
theme_bw()
limp <- read.csv("/Users/aaronsi/Documents/Documents - Aaron’s MacBook Pro/4259A/data/limpet.csv")
glimpse(limp)
## Rows: 24
## Columns: 3
## $ DENSITY <int> 8, 8, 8, 8, 8, 8, 15, 15, 15, 15, 15, 15, 30, 30, 30, 30, 30, …
## $ SEASON <chr> "spring", "spring", "spring", "summer", "summer", "summer", "s…
## $ EGGS <dbl> 2.875, 2.625, 1.750, 2.125, 1.500, 1.875, 2.600, 1.866, 2.066,…
ggplot(limp, aes(x = DENSITY, y = EGGS, colour = SEASON)) +
geom_point() +
scale_alpha_manual(values = c(spring = "green", summer = "red")) +
theme_bw()
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's alpha values.
limp.mod <- lm(EGGS ~ DENSITY * SEASON, data = limp)
names(limp.mod)
## [1] "coefficients" "residuals" "effects" "rank"
## [5] "fitted.values" "assign" "qr" "df.residual"
## [9] "contrasts" "xlevels" "call" "terms"
## [13] "model"
autoplot(limp.mod, smooth.colour = NA)
## Warning: Removed 24 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 24 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 24 rows containing missing values or values outside the scale range
## (`geom_line()`).
anova(limp.mod)
## Analysis of Variance Table
##
## Response: EGGS
## Df Sum Sq Mean Sq F value Pr(>F)
## DENSITY 1 5.0241 5.0241 30.1971 2.226e-05 ***
## SEASON 1 3.2502 3.2502 19.5350 0.0002637 ***
## DENSITY:SEASON 1 0.0118 0.0118 0.0711 0.7925333
## Residuals 20 3.3275 0.1664
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(limp.mod)
##
## Call:
## lm(formula = EGGS ~ DENSITY * SEASON, data = limp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.65468 -0.25021 -0.03318 0.28335 0.57532
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.664166 0.234118 11.380 3.45e-10 ***
## DENSITY -0.033650 0.008259 -4.074 0.000591 ***
## SEASONsummer -0.812282 0.331092 -2.453 0.023450 *
## DENSITY:SEASONsummer 0.003114 0.011680 0.267 0.792533
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4079 on 20 degrees of freedom
## Multiple R-squared: 0.7135, Adjusted R-squared: 0.6705
## F-statistic: 16.6 on 3 and 20 DF, p-value: 1.186e-05
expand.grid(FIRST = c("A", "B"), SECOND = c(1,2))
## FIRST SECOND
## 1 A 1
## 2 B 1
## 3 A 2
## 4 B 2
predict(limp.mod)
## 1 2 3 4 5 6 7 8
## 2.3949692 2.3949692 2.3949692 1.6075953 1.6075953 1.6075953 2.1594217 2.1594217
## 9 10 11 12 13 14 15 16
## 2.1594217 1.3938428 1.3938428 1.3938428 1.6546769 1.6546769 1.6546769 0.9358016
## 17 18 19 20 21 22 23 24
## 0.9358016 0.9358016 1.1499321 1.1499321 1.1499321 0.4777604 0.4777604 0.4777604
new.x <- expand.grid(DENSITY =
seq(from = 8, to = 45, length.out = 10))
head(new.x)
## DENSITY
## 1 8.00000
## 2 12.11111
## 3 16.22222
## 4 20.33333
## 5 24.44444
## 6 28.55556
limp$SEASON <- factor(limp$SEASON)
new.x <- expand.grid(
DENSITY = seq(from = 8, to = 45, length.out = 10),
SEASON = levels(limp$SEASON))
head(new.x)
## DENSITY SEASON
## 1 8.00000 spring
## 2 12.11111 spring
## 3 16.22222 spring
## 4 20.33333 spring
## 5 24.44444 spring
## 6 28.55556 spring
new.y <- predict(limp.mod, newdata = new.x,
interval = 'confidence')
head(new.y)
## fit lwr upr
## 1 2.394969 2.019285 2.770654
## 2 2.256632 1.931230 2.582034
## 3 2.118294 1.834274 2.402315
## 4 1.979957 1.724062 2.235852
## 5 1.841619 1.595998 2.087241
## 6 1.703282 1.447918 1.958646
addThese <- data.frame(new.x, new.y)
addThese <- rename(addThese, EGGS = fit)
head(addThese)
## DENSITY SEASON EGGS lwr upr
## 1 8.00000 spring 2.394969 2.019285 2.770654
## 2 12.11111 spring 2.256632 1.931230 2.582034
## 3 16.22222 spring 2.118294 1.834274 2.402315
## 4 20.33333 spring 1.979957 1.724062 2.235852
## 5 24.44444 spring 1.841619 1.595998 2.087241
## 6 28.55556 spring 1.703282 1.447918 1.958646
new.x2 <- expand.grid(DENSITY = mean(limp$DENSITY),
SEASON = levels(limp$SEASON))
predEgg <- predict(limp.mod, newdata = new.x2)
EggAtMeanDens <- data_frame(new.x2, predEgg)
## Warning: `data_frame()` was deprecated in tibble 1.1.0.
## ℹ Please use `tibble()` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
head(EggAtMeanDens)
## # A tibble: 2 × 3
## DENSITY SEASON predEgg
## <dbl> <fct> <dbl>
## 1 24.5 spring 1.84
## 2 24.5 summer 1.10
ggplot(limp, aes(x = DENSITY, y = EGGS, colour = SEASON)) +
geom_point(size = 5) +
geom_smooth(data = addThese,
aes(ymin = lwr, ymax = upr,
fill = SEASON), stat = 'identity') +
scale_colour_manual(values = c("green", summer = "red")) +
scale_fill_manual(values = c(spring = "green", summer = "red")) +
theme_bw()
