diff --git a/lessons/R-complex-time-series-models-1/material.qmd b/lessons/R-complex-time-series-models-1/material.qmd index 7a221e1..c13669d 100644 --- a/lessons/R-complex-time-series-models-1/material.qmd +++ b/lessons/R-complex-time-series-models-1/material.qmd @@ -8,7 +8,7 @@ order: 1 ### Data -Download the [Desert Pocket Mouse data](/data/pp_abundance_timeseries.csv) +Download the [Desert Pocket Mouse data](/data/pp_abundance_by_month.csv) ### Software Installation diff --git a/lessons/R-complex-time-series-models-1/r_tutorial.qmd b/lessons/R-complex-time-series-models-1/r_tutorial.qmd index 91009e3..94ddaae 100644 --- a/lessons/R-complex-time-series-models-1/r_tutorial.qmd +++ b/lessons/R-complex-time-series-models-1/r_tutorial.qmd @@ -51,7 +51,7 @@ library(dplyr) * Data on the population dynamics of the Desert Pocket Mouse ```r -pp_data <- read.csv("pp_abundance_timeseries.csv") +pp_data <- read.csv("pp_abundance_by_month.csv") ``` * mvgam doesn't currently work with tsibbles @@ -60,10 +60,12 @@ pp_data <- read.csv("pp_abundance_timeseries.csv") * Helps when analyzing multiple time series at once (e.g., multiple species) ```r -pp_data <- read.csv("pp_abundance_timeseries.csv") |> - mutate(time = newmoonnumber) |> - mutate(series = as.factor('PP')) |> - select(time, series, abundance, mintemp, cool_precip) +pp_data = read.csv("pp_abundance_by_month.csv") |> + mutate(month = yearmonth(month)) |> + mutate(time = as.numeric(month)) |> + mutate(series = as.factor('PP')) |> + tibble() +pp_data ``` * 124 months of data diff --git a/lessons/R-complex-time-series-models-2/material.qmd b/lessons/R-complex-time-series-models-2/material.qmd index becee69..f614df8 100644 --- a/lessons/R-complex-time-series-models-2/material.qmd +++ b/lessons/R-complex-time-series-models-2/material.qmd @@ -8,7 +8,7 @@ order: 1 ### Data -Download the [Desert Pocket Mouse data](/data/pp_abundance_timeseries.csv) +Download the [Desert Pocket Mouse data](/data/pp_abundance_by_month.csv) ### Software Installation diff --git a/lessons/R-complex-time-series-models-2/r_tutorial.qmd b/lessons/R-complex-time-series-models-2/r_tutorial.qmd index efaa343..b76e541 100644 --- a/lessons/R-complex-time-series-models-2/r_tutorial.qmd +++ b/lessons/R-complex-time-series-models-2/r_tutorial.qmd @@ -51,7 +51,7 @@ library(dplyr) * Data on the population dynamics of the Desert Pocket Mouse ```r -pp_data <- read.csv("pp_abundance_timeseries.csv") +pp_data <- read.csv("pp_abundance_by_month.csv") ``` * mvgam doesn't currently work with tsibbles @@ -60,10 +60,12 @@ pp_data <- read.csv("pp_abundance_timeseries.csv") * Helps when analyzing multiple time series at once (e.g., multiple species) ```r -pp_data <- read.csv("pp_abundance_timeseries.csv") |> - mutate(time = newmoonnumber) |> - mutate(series = as.factor('PP')) |> - select(time, series, abundance, mintemp, cool_precip) +pp_data = read.csv("pp_abundance_by_month.csv") |> + mutate(month = yearmonth(month)) |> + mutate(time = as.numeric(month)) |> + mutate(series = as.factor('PP')) |> + tibble() +pp_data ``` * 124 months of data @@ -91,7 +93,7 @@ $$y_t = c + \beta_1 x_{1,t} + \beta_2 y_{t-1} + \mathcal{N}(0,\sigma^{2})$$ * Which can also be written as -$$y_t = \mathcal{N}(\mu,\sigma^{2})$$ +$$y_t = \mathcal{N}(\mu_t,\sigma^{2})$$ $$u_t = c + \beta_1 x_{1,t} + \beta_2 y_{t-1}$$ * Fit using the `mvgam()` function @@ -157,14 +159,8 @@ plot(baseline_model) ### Bayesian model forecasting -* So let's look at the forecasts -* To make forecasts for Bayesian models we include data for `y` that is `NA` -* This tells the model that we don't know the values and therefore the model estimates them as part of the fitting process -* Handy because it means these models can handle missing data -* To make a true forecast one `NA` is added to the end of `y` for each time step we want to forecast -* To hindcast the values for `y` that are part of the test set are replaced with `NA` - -* We can do this automatically in mvgam using the `forecast()` function +* To make forecasts in mvgam we use the `forecast()` function (just like fable) +* Use `newdata` (not `new_data` as in fable) to include the test data for the driver forecasts ```r baseline_forecast = forecast(baseline_model, newdata = data_test) @@ -189,33 +185,23 @@ plot(baseline_forecast) $$y_t = \mathrm{Pois}(\lambda_t)$$ - * The Poisson distribution has one parameter $\lambda$ * Which is both the mean and the variance * It generates only integer draws based on a mean ```r -rpois(n = 10, lambda = 5) -hist(rpois(n = 1000, lambda = 5)) +rpois(n = 10, lambda = 4.5) +hist(rpois(n = 1000, lambda = 4.5)) ``` -* If the mean is 1.5 sometimes you'll draw a 1, sometimes a 2, sometimes a 0, etc. +* If the mean is 4.5 sometimes you'll draw a 4, sometimes a 5, sometimes a 0 or a 10 * $\lambda_t$ can be a decimal - -```r -hist(rpois(n = 1000, lambda = 4.5)) -``` - * We could expect an average of 4.5 rodents based on the environment even though we can only observe an integer number -* But $\lambda_t$ does have to be positive because we can't reasonably expect to see negative rodents - -```r -hist(rpois(n = 1000, lambda = -4.5)) -``` +* But $\lambda_t$ it does have to be positive because we can't reasonably expect to see negative rodents -* To handle this we use a log link function to give us only positive values of $\mu_t$ -* The log link means that instead of modeling $\mu_t$ directly we model $log(\mu_t)$ +* To handle this we use a log link function to give us only positive values of $\lambda_t$ +* The log link means that instead of modeling $\lambda_t$ directly we model $log(\lambda_t)$ $$y_t = \mathrm{Pois}(\lambda_t)$$ $$\mathrm{log(\lambda_t)} = c + \beta_1 x_{1,t} + \beta_2 y_{t-1}$$ @@ -326,7 +312,7 @@ plot(poisson_gam_model) * We can also see that the residual season autocorrelation is improved -## State space models +## State space models (optional) * To add an explicit observation model we use the @@ -389,5 +375,5 @@ scores <- score(poisson_gam_forecast, interval_width = 0.5) ```r in_interval = scores$PP$in_interval -length(in_interval[in_interval == TRUE]) / length(in_interval) +length(in_interval[in_interval == 1]) / length(in_interval) ``` diff --git a/lessons/R-time-series-modeling-2/material.qmd b/lessons/R-time-series-modeling-2/material.qmd index 4187dc8..b30e5c9 100644 --- a/lessons/R-time-series-modeling-2/material.qmd +++ b/lessons/R-time-series-modeling-2/material.qmd @@ -5,5 +5,7 @@ order: 1 Before starting lesson: -* Make sure the Portal time series data is available on your computer -* No new packages need to be installed for this lesson +```r +install.packages(c('tsibble', 'fable', 'feasts', 'ggtime')) +download.file("https://course.naturecast.org/data/portal_timeseries.csv", "portal_timeseries.csv") +``` diff --git a/lessons/R-time-series-modeling-3/material.qmd b/lessons/R-time-series-modeling-3/material.qmd index 5c09878..e61d9ea 100644 --- a/lessons/R-time-series-modeling-3/material.qmd +++ b/lessons/R-time-series-modeling-3/material.qmd @@ -5,6 +5,7 @@ order: 1 Before starting lesson: -* Make sure the Portal time series data is available on your computer -* Download the [Desert Pocket Mouse data](/data/pp_abundance_timeseries.csv) -* No new packages need to be installed for this lesson +```r +install.packages(c('tsibble', 'fable', 'feasts', 'ggtime')) +download.file("https://course.naturecast.org/data/pp_abundance_by_month.csv", "pp_abundance_by_month.csv") +``` diff --git a/lessons/R-time-series-modeling-3/r_tutorial.qmd b/lessons/R-time-series-modeling-3/r_tutorial.qmd index 2a7ed9b..f76cd25 100644 --- a/lessons/R-time-series-modeling-3/r_tutorial.qmd +++ b/lessons/R-time-series-modeling-3/r_tutorial.qmd @@ -21,8 +21,9 @@ library(ggtime) * Load the data ```r -pp_data = read.csv("pp_abundance_timeseries.csv") |> - as_tsibble(index = newmoonnumber) +pp_data = read.csv("pp_abundance_by_month.csv") |> + mutate(month = yearmonth(month)) |> + as_tsibble(index = month) pp_data ``` diff --git a/lessons/R-time-series-modeling/material.qmd b/lessons/R-time-series-modeling/material.qmd index 4187dc8..b30e5c9 100644 --- a/lessons/R-time-series-modeling/material.qmd +++ b/lessons/R-time-series-modeling/material.qmd @@ -5,5 +5,7 @@ order: 1 Before starting lesson: -* Make sure the Portal time series data is available on your computer -* No new packages need to be installed for this lesson +```r +install.packages(c('tsibble', 'fable', 'feasts', 'ggtime')) +download.file("https://course.naturecast.org/data/portal_timeseries.csv", "portal_timeseries.csv") +```