library(readr)
growth <- read_csv("/Users/clairedilosa/Desktop/Bio9915/Bio9915/data/datasets-master/growth.csv")
## Rows: 48 Columns: 3
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (2): supplement, diet
## dbl (1): gain
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
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)
## 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,…
levels(as.factor(growth$diet))
## [1] "barley" "oats" "wheat"
levels(as.factor(growth$supplement))
## [1] "agrimore" "control" "supergain" "supersupp"
growth <- mutate(growth, supplement = relevel(as.factor(supplement), ref="control"))
levels(as.factor(growth$supplement))
## [1] "control" "agrimore" "supergain" "supersupp"
sumMoo <- growth %>%
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
## <chr> <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)
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)
##
## 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 %>%
group_by(diet, supplement) %>%
summarize (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()+
geom_errorbar(aes(ymin=meangrow-seGrow, ymax = meangrow+seGrow, width = 0.1)) + theme_bw( )
limp <- read_csv("/Users/clairedilosa/Desktop/Bio9915/Bio9915/data/datasets-master/limpet.csv")
## Rows: 24 Columns: 3
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): SEASON
## dbl (2): DENSITY, EGGS
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
glimpse(limp)
## Rows: 24
## Columns: 3
## $ DENSITY <dbl> 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_colour_manual(values=c(spring="green", summer="red"))+ theme_bw()
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
new.x <-expand.grid(DENSITY = seq(from=8, to=45, length.out = 10), SEASON = levels(as.factor(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.x <-expand.grid(DENSITY=mean(limp$DENSITY), SEASON=levels(as.factor(limp$SEASON)))
predEgg <-predict(limp.mod,newdata = new.x)
EggAtMeanDens <- data.frame(new.x, predEgg)
head(EggAtMeanDens)
## DENSITY SEASON predEgg
## 1 24.5 spring 1.83975
## 2 24.5 summer 1.10375
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(spring='green', summer='red'))+ scale_fill_manual(values=c(spring='green', summer='red'))+theme_bw()
