Project 2
Natalia Sifnugel
2023-09-24
library(tidyverse)
library(tidytuesdayR)
library(here)
library(lubridate)
library(sessioninfo)Part 1 Fun with functions
Part 1A: Exponential transformation
Exp <- function(x, k) {
exp = 1
for(i in 1:k){
exp2 = x^i/factorial(i)
exp = exp + exp2
}
return(exp)
}Sample outputs
Exp(2,4)## [1] 7# check
1 + 2 + (2^2/factorial(2)) + (2^3/factorial(3)) + (2^4/factorial(4))## [1] 7Part 1B: Sample mean and sample standard deviation
sample_mean <- function(x) {
numerator = 0
denominator = 0
for(i in seq_along(x)){
numerator = numerator + x[i]
denominator = denominator + 1
}
mean = numerator/denominator
return(mean)
}
sample_sd <- function(x) {
numerator = 0
denominator = 0
for(i in seq_along(x)){
numerator = numerator + (x[i]-sample_mean(x))^2
denominator = denominator + 1
}
sd = sqrt(numerator/(denominator-1))
return(sd)
}Sample outputs
x = c(1,22,14,67,10)
sample_mean(x)## [1] 22.8# check
mean(x)## [1] 22.8sample_sd(x)## [1] 25.83989# check
sd(x)## [1] 25.83989Part 1C: Confidence intervals
calculate_CI <- function(x, conf = 0.95) {
alpha <- 1 - conf
degrees_freedom <- length(x) - 1
t_score <- qt(p = alpha / 2, df = degrees_freedom, lower.tail = FALSE)
lower_tail = sample_mean(x) - t_score*(sample_sd(x)/sqrt(length(x)))
upper_tail = sample_mean(x) + t_score*(sample_sd(x)/sqrt(length(x)))
CI = c(lower_tail,upper_tail)
return(CI)
}Sample outputs
x = c(1,22,14,67,10)
calculate_CI(x,0.90)## [1] -1.835517 47.435517calculate_CI(x,0.80)## [1] 5.082344 40.517656# check
dat <- data.frame(x = x)
fit <- lm(x ~ 1, dat)
confint(fit, level = 0.90)## 5 % 95 %
## (Intercept) -1.835517 47.43552confint(fit, level = 0.80)## 10 % 90 %
## (Intercept) 5.082344 40.51766Part 2: Wrangling data
Downloading datasets
if (!dir.exists(here("data"))) {
dir.create(here("data"))
}
if (!file.exists(here("data", "tuesdata_rainfall.RDS"))) {
tuesdata <- tidytuesdayR::tt_load("2020-01-07")
rainfall <- tuesdata$rainfall
temperature <- tuesdata$temperature
# save the files to RDS objects
saveRDS(tuesdata$rainfall, file = here("data", "tuesdata_rainfall.RDS"))
saveRDS(tuesdata$temperature, file = here("data", "tuesdata_temperature.RDS"))
}
rainfall <- readRDS(here("data", "tuesdata_rainfall.RDS"))
temperature <- readRDS(here("data", "tuesdata_temperature.RDS"))Wrangling rainfall dataset and joining it to temperature dataset
rainfall_temp <-rainfall %>%
drop_na() %>%
mutate(date = ymd(paste(year,month,day, sep = "-")), city_name = toupper(city_name)) %>%
select(-month,-day) %>%
inner_join(temperature, by = c("date","city_name"), relationship = "many-to-many")
show(rainfall_temp)## # A tibble: 83,964 × 13
## station_code city_name year rainfall period quality lat long station_name
## <chr> <chr> <dbl> <dbl> <dbl> <chr> <dbl> <dbl> <chr>
## 1 009151 PERTH 1967 2.8 1 Y -32.0 116. Subiaco Was…
## 2 009151 PERTH 1967 2.8 1 Y -32.0 116. Subiaco Was…
## 3 009151 PERTH 1967 4.8 1 Y -32.0 116. Subiaco Was…
## 4 009151 PERTH 1967 4.8 1 Y -32.0 116. Subiaco Was…
## 5 009151 PERTH 1967 5.8 1 Y -32.0 116. Subiaco Was…
## 6 009151 PERTH 1967 5.8 1 Y -32.0 116. Subiaco Was…
## 7 009151 PERTH 1967 16 1 Y -32.0 116. Subiaco Was…
## 8 009151 PERTH 1967 16 1 Y -32.0 116. Subiaco Was…
## 9 009151 PERTH 1967 1 1 Y -32.0 116. Subiaco Was…
## 10 009151 PERTH 1967 1 1 Y -32.0 116. Subiaco Was…
## # ℹ 83,954 more rows
## # ℹ 4 more variables: date <date>, temperature <dbl>, temp_type <chr>,
## # site_name <chr>Part 3: Data visualization
Part 3A: Plotting temperature data over time
rainfall_temp %>%
filter(year >= 2014) %>%
ggplot() +
geom_line(aes(x = date, y = temperature, color = temp_type))+
facet_grid(city_name ~. ) +
labs(title = "Temperature vs. Year",
subtitle = "Canberra experiences greatest range in daily temperature. In Sydney, peak daily temperatures are rising.",
caption = "Source: Australian Government, Bureau of Meteorology ",
x = "Date",
y = "Temperature (C)") +
theme_light() +
theme(axis.title.x = element_text(size = 16),
axis.title.y = element_text(size = 16),
plot.title = element_text(size = 22),
plot.subtitle = element_text(size = 16))
Part 3B: Plotting rainfall over time
rainfall_histogram <- function(x, y) {
if(x %in% rainfall_temp$city_name){
rainfall_city <- rainfall_temp %>%
filter(city_name == x)
} else {
stop(x," is not in the dataset.")
}
if(y %in% rainfall_city$year){
rainfall_city %>%
filter(year == y) %>%
ggplot(aes(log(rainfall + 0.000000000000001))) +
geom_histogram(fill = 'lightblue', color = 'black') +
labs(title = paste("Distribution of Rainfall in",x,"in",y, sep = " "),
subtitle = paste("Use this histogram to determine if",y,"was a year of deluge, drought, or expected rain in",x,sep = " "),
caption = "Source: Australian Government, Bureau of Meteorology ",
x = "log(Rainfall in mm)",
y = "Number of Days") +
theme_light() +
theme(axis.title.x = element_text(size = 12),
axis.title.y = element_text(size = 12),
plot.title = element_text(size = 18),
plot.subtitle = element_text(size = 12))
} else{
stop(y," not in dataset for ",x,".")
}
}Function design
In this function, I used an if statement to check if the city is in the data set. If the city is in the data set, then a new filtered data frame is created in which only the values for the specified city are included. This data set then enters a new if else statement to check if the year (second input) exists for this city. If the year does exist in the filtered city data set, then a histogram is generated for that city for that year. If these two conditions were not met, the function returns a statement indicating the data are not present.
Because the log(0) is an infinite number, days with 0 mm of rainfall were initially dropped from the histogram. Days of 0 rainfall for this data set as a whole are important because it seeks to understand the relationship between temperate, rain, and fires. So, to ensure that days of 0 rainfall are included, I added a small constant (rainfall + 0.000000000000001).
Example output
rainfall_histogram('PERTH',2000)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
rainfall_histogram('BRISBANE',1981)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
rainfall_histogram('SYDNEY',2014)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Part 4: Apply functions and plot
Part 4A: Summary of rainfall data for each city since 2014
rain_df <- rainfall_temp %>%
filter(year >= 2014) %>%
group_by(city_name, year) %>%
summarize(tibble(
`mean` = sample_mean(rainfall),
`sd` = sample_sd(rainfall),
`lower_bound` = calculate_CI(rainfall,0.95)[1],
`upper_bound` = calculate_CI(rainfall,0.95)[2]
)
)## `summarise()` has grouped output by 'city_name'. You can override using the
## `.groups` argument.show(rain_df)## # A tibble: 30 × 6
## # Groups: city_name [5]
## city_name year mean sd lower_bound upper_bound
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 BRISBANE 2014 3.84 12.3 2.99 4.69
## 2 BRISBANE 2015 6.17 22.7 4.53 7.82
## 3 BRISBANE 2016 2.09 8.61 1.47 2.72
## 4 BRISBANE 2017 2.72 10.2 1.97 3.47
## 5 BRISBANE 2018 2.47 9.63 1.75 3.19
## 6 BRISBANE 2019 2.10 6.35 1.38 2.82
## 7 CANBERRA 2014 1.57 4.84 1.22 1.92
## 8 CANBERRA 2015 1.51 5.16 1.13 1.89
## 9 CANBERRA 2016 2.15 6.28 1.70 2.61
## 10 CANBERRA 2017 1.34 4.12 1.04 1.64
## # ℹ 20 more rowsPart 4B: Plotting mean rainfall
rain_df %>%
ggplot() +
geom_point(aes(year,mean, color = city_name)) +
geom_line(aes(year,mean, color = city_name)) +
facet_grid(.~ city_name) +
geom_errorbar(aes(ymin= lower_bound, ymax= upper_bound, x = year), width = 0.2) +
labs(title = "Mean Rainfall by Year",
subtitle = "Average rainfall is on the decline across major cities in Australia.",
caption = "Source: Australian Government, Bureau of Meteorology ",
x = "Year",
y = "Rainfall in mm") +
theme_light() +
theme(legend.position = "none",
axis.title.x = element_text(size = 16),
axis.title.y = element_text(size = 16),
plot.title = element_text(size = 22),
plot.subtitle = element_text(size = 16))
R session information
options(width = 120)
sessioninfo::session_info()## ─ Session info ───────────────────────────────────────────────────────────────────────────────────────────────────────
## setting value
## version R version 4.3.1 (2023-06-16)
## os macOS Ventura 13.4.1
## system aarch64, darwin20
## ui X11
## language (EN)
## collate en_US.UTF-8
## ctype en_US.UTF-8
## tz America/New_York
## date 2023-09-27
## pandoc 3.1.1 @ /Applications/RStudio.app/Contents/Resources/app/quarto/bin/tools/ (via rmarkdown)
##
## ─ Packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────
## package * version date (UTC) lib source
## bslib 0.5.1 2023-08-11 [1] CRAN (R 4.3.0)
## cachem 1.0.8 2023-05-01 [1] CRAN (R 4.3.0)
## cellranger 1.1.0 2016-07-27 [1] CRAN (R 4.3.0)
## cli 3.6.1 2023-03-23 [1] CRAN (R 4.3.0)
## colorspace 2.1-0 2023-01-23 [1] CRAN (R 4.3.0)
## curl 5.0.2 2023-08-14 [1] CRAN (R 4.3.0)
## digest 0.6.33 2023-07-07 [1] CRAN (R 4.3.0)
## dplyr * 1.1.3 2023-09-03 [1] CRAN (R 4.3.0)
## evaluate 0.21 2023-05-05 [1] CRAN (R 4.3.0)
## fansi 1.0.4 2023-01-22 [1] CRAN (R 4.3.0)
## farver 2.1.1 2022-07-06 [1] CRAN (R 4.3.0)
## fastmap 1.1.1 2023-02-24 [1] CRAN (R 4.3.0)
## forcats * 1.0.0 2023-01-29 [1] CRAN (R 4.3.0)
## fs 1.6.3 2023-07-20 [1] CRAN (R 4.3.0)
## generics 0.1.3 2022-07-05 [1] CRAN (R 4.3.0)
## ggplot2 * 3.4.3 2023-08-14 [1] CRAN (R 4.3.0)
## glue 1.6.2 2022-02-24 [1] CRAN (R 4.3.0)
## gtable 0.3.4 2023-08-21 [1] CRAN (R 4.3.0)
## here * 1.0.1 2020-12-13 [1] CRAN (R 4.3.0)
## hms 1.1.3 2023-03-21 [1] CRAN (R 4.3.0)
## htmltools 0.5.6 2023-08-10 [1] CRAN (R 4.3.0)
## httr 1.4.7 2023-08-15 [1] CRAN (R 4.3.0)
## jquerylib 0.1.4 2021-04-26 [1] CRAN (R 4.3.0)
## jsonlite 1.8.7 2023-06-29 [1] CRAN (R 4.3.0)
## knitr 1.44 2023-09-11 [1] CRAN (R 4.3.0)
## labeling 0.4.3 2023-08-29 [1] CRAN (R 4.3.0)
## lifecycle 1.0.3 2022-10-07 [1] CRAN (R 4.3.0)
## lubridate * 1.9.2 2023-02-10 [1] CRAN (R 4.3.0)
## magrittr 2.0.3 2022-03-30 [1] CRAN (R 4.3.0)
## munsell 0.5.0 2018-06-12 [1] CRAN (R 4.3.0)
## pillar 1.9.0 2023-03-22 [1] CRAN (R 4.3.0)
## pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.3.0)
## purrr * 1.0.2 2023-08-10 [1] CRAN (R 4.3.0)
## R6 2.5.1 2021-08-19 [1] CRAN (R 4.3.0)
## readr * 2.1.4 2023-02-10 [1] CRAN (R 4.3.0)
## readxl 1.4.3 2023-07-06 [1] CRAN (R 4.3.0)
## rlang 1.1.1 2023-04-28 [1] CRAN (R 4.3.0)
## rmarkdown 2.24 2023-08-14 [1] CRAN (R 4.3.0)
## rprojroot 2.0.3 2022-04-02 [1] CRAN (R 4.3.0)
## rstudioapi 0.15.0 2023-07-07 [1] CRAN (R 4.3.0)
## rvest 1.0.3 2022-08-19 [1] CRAN (R 4.3.0)
## sass 0.4.7 2023-07-15 [1] CRAN (R 4.3.0)
## scales 1.2.1 2022-08-20 [1] CRAN (R 4.3.0)
## sessioninfo * 1.2.2 2021-12-06 [1] CRAN (R 4.3.0)
## stringi 1.7.12 2023-01-11 [1] CRAN (R 4.3.0)
## stringr * 1.5.0 2022-12-02 [1] CRAN (R 4.3.0)
## tibble * 3.2.1 2023-03-20 [1] CRAN (R 4.3.0)
## tidyr * 1.3.0 2023-01-24 [1] CRAN (R 4.3.0)
## tidyselect 1.2.0 2022-10-10 [1] CRAN (R 4.3.0)
## tidytuesdayR * 1.0.2 2022-02-01 [1] CRAN (R 4.3.0)
## tidyverse * 2.0.0 2023-02-22 [1] CRAN (R 4.3.0)
## timechange 0.2.0 2023-01-11 [1] CRAN (R 4.3.0)
## tzdb 0.4.0 2023-05-12 [1] CRAN (R 4.3.0)
## usethis 2.2.2 2023-07-06 [1] CRAN (R 4.3.0)
## utf8 1.2.3 2023-01-31 [1] CRAN (R 4.3.0)
## vctrs 0.6.3 2023-06-14 [1] CRAN (R 4.3.0)
## withr 2.5.0 2022-03-03 [1] CRAN (R 4.3.0)
## xfun 0.40 2023-08-09 [1] CRAN (R 4.3.0)
## xml2 1.3.5 2023-07-06 [1] CRAN (R 4.3.0)
## yaml 2.3.7 2023-01-23 [1] CRAN (R 4.3.0)
##
## [1] /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library
##
## ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────