19.3 — Monte Carlo methods

edits: — under construction —

Introduction

Statistical method that employ Monte Carlo methods use repeated random sampling to estimate properties of a frequency distribution. These distributions may be well-known, e.g., gamma-distribution, normal distribution, or t-distribution, or  . The simulation is based on generation of a set of random numbers on the open interval (0,1) — the set of real numbers between zero and one (all numbers greater than 0 and less than 1).

Markov chain Monte Carlo (MCMC) sampling approach used to solve large scale problems. The Markov chain refers to how the sample is drawn from a specified probability distribution. It can be drawn by discrete time steps (DTMC) or by a continuous process (CTMC). The Markov process is “memoryless:” predictions of future events are derived solely from their present state — the future and past states are independent.

Gibbs sampling is a common MCMC algorithm.

ccc

R code

R’s uniform generator is runif function. Examples of the samples generated over different values (100, 1000, 10000, 100000) with output displayed as histograms (Fig. 1). Note that as sample size increases, the simulated distributions resemble more and more the uniform distribution. Use set.seed() to reproduce the same set and sequence of numbers

require(RcmdrMisc)
par(mfrow = c(2, 2))
myUniformH <- data.frame(runif(100))
with(myUniformH, Hist(runif.100., scale="frequency", ylim=c(0,20), breaks="Sturges", col="red", xlab="100 samples", ylab="Count"))
myUniform1K <- data.frame(runif(1000))
with(myUniform1K, Hist(runif.1000., scale="frequency", ylim=c(0,150), breaks="Sturges", col="green", xlab="1K samples", ylab="Count"))
myUniform10K <- data.frame(runif(10000))
with(myUniform10K, Hist(runif.10000., scale="frequency", ylim=c(0,600), breaks="Sturges", col="lightblue", xlab="10K samples", ylab="Count"))
myUniform100K <- data.frame(runif(100000))
with(myUniform100K, Hist(runif.100000., scale="frequency", ylim=c(0,5000),breaks="Sturges", col="blue", xlab="100K samples", ylab="Count"))
#reset par()
dev.off()

Figure 1. Histograms of runif results with 100, 1K, 10K, and 100K numbers of values to be generated 

Looks pretty uniform. A property of random numbers is that history should not influence the future, i.e., no autocorrelation. We can check using the acf() function (Fig. 2).

par(mfrow = c(2, 2)) 
acf(myUniformH, main="100")
acf(myUniform1K, main="1K")
acf(myUniform10K, main="10K")
acf(myUniform100K, main="100K"
dev.off()

Figure 2. Autocorrelation plots of runif results with 100, 1K, 10K, and 100K numbers of values

Correlations among points are plotted versus lag, where lag refers to the number of points between adjacent points, e.g., lag = 10 reflects the correlation among points 1 and 11, 2 and 12, and so forth. The band defined by two parallel blue dashed lines

Questions

1. Use set.seed(123) and repeat runif(10) twice. Confirm that the two sets are different (do not set seed) or the same when set.seed is used. R hint: use function identical(x,y), where x and y are the two generated samples. This function tests whether the values and sequence of elements are the same between the two vectors.

Chapter 19 contents

• Introduction
• Jackknife sampling
• Bootstrap sampling
• Monte Carlo methods
• Ch19 References and suggested readings

17.8 – Assumptions and model diagnostics for Simple Linear Regression

Introduction

The assumptions for all linear regression:

1. Linear model is appropriate.
   The data are well described (fit) by a linear model.
2. Independent values of Y and equal variances.
   Although there can be more than one Y for any value of X, the Y‘s cannot be related to each other (that’s what we mean by independent). Since we allow for multiple Y‘s for each X, then we assume that the variances of the range of Y‘s are equal for each X value (this is similar to our ANOVA assumptions for equal variance by groups). Another term for equal variances is homoscedasticity.
3. Normality.
   For each X value there is a normal distribution of Y‘s (think of doing the experiment over and over).
4. Error
   The residuals (error) are normally distributed with a mean of zero.

Note the mnemonic device: Linear, Independent, Normal, Error or LINE.

Each of the four elements will be discussed below in the context of Model Diagnostics. These assumptions apply to how the model fits the data. There are other assumptions that, if violated, imply you should use a different method for estimating the parameters of the model.

Ordinary least squares makes the additional assumption about the quality of the independent variable that e that measurement of X is done without error. Measurement error is a fact of life in science, but the influence of error on regression differs if the error is associated with the dependent or independent variable. Measurement error in the dependent variable increases the dispersion of the residuals but will not affect the estimates of the coefficients; error associated with the independent variables, however, will affect estimates of the slope. In short, error in X leads to biased estimates of the slope.

The equivalent, but less restrictive practical application of this assumption is that the error in X is at least negligible compared to the measurements in the dependent variable.

Multiple regression makes one more assumption, about the relationship between the predictor variables (the X variables). The assumption is that there is no multicollinearity, a subject we will bring up next time (see Chapter 18).

Model diagnostics

We just reviewed how to evaluate the estimates of the coefficients of the model. Now we need to address a deeper meaning — how well the model explains the data. Consider a simple linear regression first. If H_o: b = 0 is not rejected — then the slope of the regression equation is taken to not differ from zero. We would conclude that if repeated samples were drawn from the population, on average, the regression equation would not fit the data well (lots of scatter) and it would not yield useful prediction.

However, recall that we assume that the fit is linear. One assumption we make in regression is that a line can, in fact, be used to describe the relationship between X and Y.

Here are two very different situations where the slope = 0.

Example 1. Linear Slope = 0, No relationship between X and Y

Example 2. Linear Slope = 0, A significant relationship between X and Y

But even if H_o: b = 0 is rejected (and we conclude that a linear relationship between X and Y is present), we still need to be concerned about the fit of the line to the data — the relationship may be more nonlinear than linear, for example. Here are two very different situations where the slope is not equal to 0.

Example 3. Linear Slope > 0, a linear relationship between X and Y

Example 4. Linear Slope > 0, curve-linear relationship between X and Y

How can you tell the difference? There are many regression diagnostic tests, many more than we can cover, but you can start with looking at the coefficient of determination (low R² means low fit to the line), and we can look at the pattern of residuals plotted against the either the predicted values or the X variables (my favorite). The important points are:

1. In linear regression, you fit a model (the slope + intercept) to the data;
2. We want the usual hypothesis tests (are the coefficients different from zero?) and
3. We need to check to see if the model fits the data well. Just like in our discussions of chi-square, a “perfect fit would mean that the difference between our model and the data would be zero.

Graph options

Using residual plots to diagnose regression equations

Yes, we need to test the coefficients (intercept H_o = 0; slope H_o = 0) of a regression equation, but we also must decide if a regression is an appropriate description of the data. This topic includes the use of diagnostic tests in regression. We address this question chiefly by looking at

1. scatterplots of the independent (predictor) variable(s) vs. dependent (response) variable(s).
   what patterns appear between X and Y? Do your eyes tell you “Line”? “Curve”? “No relation”?
2. coefficient of determination
   closer to zero than to one?
3. patterns of residuals plotted against the X variables (other types of residual plots are used to, this is one of my favorites)

Our approach is to utilize graphics along with statistical tests designed to address the assumptions.

One typical choice is to see if there are patterns in the residual values plotted against the predictor variable. If the LINE assumptions hold for your data set, then the residuals should have a mean of zero with scatter about the mean. Deviations from LINE assumptions will show up in residual plots.

Here are examples of POSSIBLE outcomes:

Figure 1. An ideal plot of residuals

Solution: Proceed! Assumptions of linear regression met.

Compare to plots of residuals that differ from the ideal.

Figure 2. We have a problem. Residual plot shows a “funnel” shape — unequal variance (aka heteroscedasticity).

Solution. Try a transform like the log₁₀-transform.

Figure 3. Problem. Residual plot shows systematic trend.

Solution. Linear model a poor fit; May be related to measurement errors for one or more predictor variables. Try adding an additional predictor variable or model the error in your general linear model.

Figure 4. Problem. Residual plot shows nonlinear trend.

Solution. Transform data or use more complex model

This is a good time to mention that in statistical analyses, one often needs to do multiple rounds of analyses, involving description and plots, tests of assumptions, tests of inference. With regression, in particular, we also need to decide if our model (e.g., linear equation) is a good description of the data.

Diagnostic plot examples

Return to our example.Tadpole dataset. To obtain residual plots, Rcmdr: Models → Graphs → Basic diagnostic plots yields four graphs. Graph A is

Figure 5. Basic diagnostic plots. A: residual plot; B: Q-Q plot of residuals; C: Scale-location (aka spread-location) plot; D: leverage residual plot.

In brief, we look at plots

A, the residual plot, to see if there are trends in the residuals. We are looking for a spread of points equally above and below the mean of zero. In Figure 5 we count seven points above and six points below zero so there’s no indication of a trend in the residuals vs the fitted VO2 (Y) values.

B, the Q-Q plot is used to see if normality holds. As discussed before, if our data are more or less normally distributed, then points will fall along a straight line in a Q-Q plot.

C, the Scale- or spread-location plot is used to verify equal variances of errors.

D, Leverage plot — looks to see if an outlier has leverage on the fit of the line to the data, i.e., changes the slope. Additionally, provides location of Cook’s distance measure (dashed red lines). Cook’s distance measures the effect on the regression by removing one point at a time and then fitting a line to the data. Points outside the dashed lines have influence.

Statistical test options

After building linear models, run statistical diagnostic tests that compliment graphics approaches. Here, simply to introduce the tests, I included the Gosner factor as potential predictor — many of these numerical diagnostics are more appropriate for models with more than one predictor variable, the subject of the next chapter.

Results for the model were:

lm(formula = VO2 ~ Body.mass + Gosner, data = frogs)

Residuals:
    Min      1Q Median    3Q    Max 
-163.12 -125.53 -20.27 83.71 228.56

Coefficients:
           Estimate Std. Error t value  Pr(>|t|) 
(Intercept) -595.37     239.87  -2.482 0.03487 * 
Body.mass    431.20     115.15   3.745 0.00459 **
Gosner[T.II]  64.96     132.83   0.489 0.63648 
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 153.4 on 9 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared: 0.8047, Adjusted R-squared: 0.7613 
F-statistic: 18.54 on 2 and 9 DF, p-value: 0.0006432

Question. Note the p-value for the coefficients — do we include a model with both body mass and Gosner Stage? See Chapter 18.5 – Selecting the best model

In R Commander, the diagnostic test options are available via

Rcmdr: Models → Numerical diagnostics

Variance inflation factors (VIF):  used to detect multicollinearity among the predictor variables. If correlations are present among the predictor variables, then you can’t rely on the the coefficient estimates — whether predictor A causes change in the response variable depends on whether the correlated B predictor is also included in the model. If correlation between predictor A and B, the statistical effect is increased variance associated with the error of the coefficient estimates. There are VIF for each predictor variable. A VIF of one means there is no correlation between that predictor and the other predictor variables. A VIF of 10 is taken as evidence of serious multicollinearity in the model.

R code:

vif(RegModel.2)

R output

Body.mass Gosner 
2.186347 2.186347 

round(cov2cor(vcov(RegModel.2)), 3) # Correlations of parameter estimates
         (Intercept) Body.mass Gosner[T.II]
(Intercept)    1.000    -0.958        0.558
Body.mass     -0.958     1.000       -0.737
Gosner[T.II]   0.558    -0.737        1.000

Our interpretation

The output shows a matrix of correlations; the correlation between the two predictor variables in this model was 0.558, large per our discussions in Chapter 16.1 . Not really surprising as older tadpoles (Gosner stage II), are larger tadpoles. We’ll leave the issue of multicollinearity to Chapter 18.1.

Breusch-Pagan test for heteroscedasticity… Recall that heteroscedasticity is another name for unequal variances. This tests pertains to the residuals; we used the scale (spread) location plot (C) above. The test statistic can be calculated as . The function call is part of the car package, which we download as part of the Rcmdr package. R Commander will pop up a menu of choices for this test. For starters, accept the defaults, enter nothing, and click OK.

R code:

library(zoo, pos=16)
library(lmtest, pos=16)
bptest(VO2 ~ Body.mass + Gosner, varformula = ~ fitted.values(RegModel.2), studentize=FALSE, data=frogs)

R output

Breusch-Pagan test

data:  VO2 ~ Body.mass + Gosner
BP = 0.0040766, df = 1, p-value = 0.9491

Our interpretation

A p-value less than 5% would be interpreted as  statistically significant unequal variances among the residuals.

Durbin-Watson for autocorrelation… used to detect serial (auto)correlation among the residuals, specifically, adjacent residuals. Substantial autocorrelation can mean Type I error for one or more model coefficients. The test returns a test statistic that ranges from 0 to 4. A value at midpoint (2) suggests no autocorrelation. Default test for autocorrelation is r > 0.

R code:

dwtest(VO2 ~ Body.mass + Gosner, alternative="greater", data=frogs)

R output

Durbin-Watson test

data:  VO2 ~ Body.mass + Gosner
DW = 2.1336, p-value = 0.4657
alternative hypothesis: true autocorrelation is greater than 0

Our interpretation

The DW statistic was about 2. Per our definition, we conclude no correlations among the residuals.

RESET test for nonlinearity…  The Ramsey Regression Equation Specification Error Test is used to interpret whether or not the equation is properly specified. The idea is that fit of the equation is compared against alternative equations that account for possible nonlinear dependencies among predictor and dependent variables. The test statistic follows an F distribution.

R code:

resettest(VO2 ~ Body.mass + Gosner, power=2:3, type="regressor", data=frogs)

R output

RESET test

data:  VO2 ~ Body.mass + Gosner
RESET = 0.93401, df1 = 2, df2 = 7, p-value = 0.437

Our interpretation

Also a nonsignificant test — adding a power term to account for a nonlinear trend gained no model fit benefits, p-value much greater than 5%.

Bonferroni outlier test. Reports the Bonferroni p-values for testing each observation to see if it can be considered to be a mean-shift outlier (a method for identifying potential outlier values — each value is in turn replaced with the mean of its nearest neighbors, Lehhman et al 2020).  Compare to results of VIF. The objective of outlier and VIF tests is to verify that the regression model is stable, insensitive to potential leverage (outlier) values.

R code:

outlierTest(LinearModel.3)

R output

No Studentized residuals with Bonferroni p < 0.05
Largest |rstudent|:
  rstudent unadjusted p-value Bonferroni p
7 1.958858        0.085809              NA

Our interpretation

We ran the outlierTest() function on LinearModel.3 to check for influential outliers using Studentized residuals with Bonferroni correction.

The output shows “No Studentized residuals with Bonferroni p < 0.05”, which means no data points were statistically flagged as outliers after adjusting for multiple comparisons. The largest Studentized residual was 1.96 (point 7), with an unadjusted p-value of 0.086. Although this is the largest deviation from the regression line, it is not statistically significant. Therefore, we conclude that no single observation is exerting undue influence on the model fit, and we can be confident that our regression results are not driven by outliers.

Response transformation… The concept of transformations of data to improve adherence to statistical assumptions was introduced in Chapters 3.2, 13.1, and 15.

R code:

summary(powerTransform(LinearModel.3, family="bcPower"))

R output

bcPower Transformation to Normality 
    Est Power Rounded Pwr Wald Lwr Bnd Wald Upr Bnd
Y1     0.9987           1       0.2755       1.7218

Likelihood ratio test that transformation parameter is equal to 0
   (log transformation)
                           LRT df      pval
LR test, lambda = (0) 7.374316  1 0.0066162

Likelihood ratio test that no transformation is needed
                                LRT df    pval
LR test, lambda = (1) 0.00001289284  1 0.99714

Our interpretation

The Box-Cox analysis suggests that the response variable for LinearModel.3 can remain on its original scale. No power transformation was needed, and any standard linear regression inference is appropriate. The estimated transformation parameter was close to one, which suggests that the original response variable did not need a transformation. The LRT tests whether a log transformation (lambda = 0) would improve normality. The significant p-value, 0.0066, for the LRT indicates that a log transformation would not be the best choice, consistent with lambda ≈ 1. And last, the very high p-value, 0.99714, confirms that no transformation is necessary; the model residuals were close to normal.

Questions

1. Write up three learning outcomes for this page. Hint: Point your favorite generative AI to this page and ask for help.
2. Referring to Figures 1 – 4 on this page, which plot best suggests a regression line fits the data?
3. Return to the electoral college data set and your linear models of Electoral vs. POP_2010 and POP_2019. Obtain the four basic diagnostic plots and comment on the fit of the regression line to the electoral college data.
   • residual plot
   • Q-Q plot
   • Scale-location plot
   • Leverage plot
4. With respect to your answers in question 2, how well does the electoral college system reflect the principle of one person one vote?

Quiz Chapter 17.8

Assumptions and model diagnostics for Simple Linear Regression

Start Quiz

Question

Confirm

✕

✓

Your answer:

Correct answer:

Next

Please wait..

↺

You got {{SCORE_CORRECT}} out of {{SCORE_TOTAL}}

Your Answers

Chapter 17 contents

• Introduction
• Simple Linear Regression
• Relationship between the slope and the correlation
• Estimation of linear regression coefficients
• OLS, RMA, and smoothing functions
• Testing regression coefficients
• Regression model fit
• Assumptions and model diagnostics for Simple Linear Regression
• References and suggested readings (Ch17 & 18)

20.5 – Time series

Introduction

Time series refers to any measure recorded over time. Stationary time series do not have trends or seasonality, just random (white) noise; differencing, or nonstationary time series do have trends and or seasonality. Stationary time series will not have predictable patterns over the long term, but for us, what this distinction says is that statistics like autocorrelation, the mean and standard deviation remains the same over time. Nonstationary times series can be made stationary by transformation, with taking a differencing approach, subtracting each observation from its preceding one to remove trends and seasonality, one common approach.

All kinds of examples of time series analysis can be made, with stock market price fluctuations (Hamilton 1994) and meteorological forecasting (Mudelsee 2014) as common examples. Time series analysis are also important in clinical situations, for example, analysis of temporal patterns in ECG and blood pressure data, which helps predict risks, monitor patient health, and understand disease trends. In ecology, perhaps the most famous time series data set is the Hudson’s Bay Company fur-trade records, which show a 9- to 11-year predator-prey population cycle between populations of the Canadian lynx (the predator) and the snowshoe hare (the prey) (Krebs et al 1995).

Significant autocorrelation — the correlation between any two observations in a data series separated by a specific time interval or lag — represents the extent to which a time series is correlated with its past values. For example, autocorrelation could measure the association between a person’s heart rate today and its heart rate yesterday, or between its heart rate today and the heart rate from last week. It is not simply the relationship between today’s heart rate and, for instance, today’s environmental temperature. This kind of comparison, between two sets of time series, calls for a cross-correlation or predictive correlation. A cross-correlation measures the correlation between the first series (temperature) and a lagged version of the second series (heart rate some minutes or hours later). Examples of such correlations include our familiar product moment correlation, but generalized across all possible time lags.

There’s much to the analysis of time series, but one key concept is the moving average. Along with any trend or pattern due to seasonality, time series data are expected to exhibit noise or random variation associated with each point in the series. The moving average is an attempt to smooth out data to reveal underlying trends and are fundamental to seasonal decomposition and forecasting. We calculate a series of averages of different subsets of the data, making it easier to spot long-term patterns by filtering out short-term noise, like fluctuations or outliers.

For example, to calculate the 5-beat moving average of a heartbeat, you would average the last 5 individual heartbeat values. As each new heartbeat is recorded, the oldest one is dropped, and a new average is calculated, creating a smoother trend line that shows a clearer overall picture rather than focusing on individual, potentially spiky, beat-to-beat variations. Rate-responsive or adaptive pacemakers use moving averages to determine the appropriate pacing rate. This helps provide a smoother and more stable heart rate response to the patient’s activity.

This statistical approach should sound familiar — we introduced without much fanfare use of LOESS (Locally Estimated Scatterplot Smoothing) to smooth data to improve pattern detection in Chapter 4.5 and Chapter 17.4. While a moving average takes average of a fixed number of recent data points, LOESS is a more powerful non-parametric regression technique that fits local polynomial regressions to subsets of the data to create a smooth curve.

In time series analysis, the autocorrelation function (ACF) measures the linear relationship between a time series and a lagged version of itself. It quantifies how correlated a series is with its past values () and ), where is the lag or time interval. An ACF plot, or correlogram, visually displays these correlation coefficients at different lags to help identify patterns, assess randomness, and determine appropriate time series models like ARIMA. ACF is related to moving averages because the ACF plot helps to identify the order of a moving average model, which depends on the ACF’s behavior at different lags.

Like many pages in Mike’s Biostatistics Book, this page only begins to introduce the subject. For a more thorough introduction, see Introduction to Time Series Analysis, NIST. See also the text by Shumway and Stoffer (2025).

ARIMA and ARMA models

ARIMA and ARMA models are fundamental tools for analyzing time series data in biostatistics, especially when the goal is to understand temporal patterns or forecast future values. An ARMA model combines two ideas: autoregression (AR), where the current value depends on previous values, and moving average (MA), where the current value depends on previous error terms. An ARIMA model extends this by adding differencing—the “I” for Integrated—which helps remove trends and make the series stationary. These models are widely used in public health, epidemiology, and biological monitoring because many variables (e.g., infection rates, physiological signals) show correlation over time.

Both ARMA and ARIMA models rely on the concepts of lagged dependence and error smoothing, but ARIMA is more flexible because it explicitly accommodates nonstationary data. An ARMA(p, q) model is appropriate when the time series is already stationary or has been preprocessed to remove trend and seasonality. is the autoregressive order, the number of lagged observations used to forecast the current value, and is the moving average order, the number of lagged forecast errors included in the model.

In contrast, ARIMA(p, d, q) includes a differencing parameter that automatically transforms the series to stationarity before fitting the ARMA structure. Thus, ARIMA models are typically used for data with a clear trend or evolving mean, while ARMA models are a subset applied to stable series without differencing. In practice, ARIMA is more commonly applied because real-world biological and clinical time series often exhibit trends.

R code

In R, ARMA and ARIMA models can be fit using the forecast or stats packages. The code below shows simple examples of each using a generic time series object . An ARMA(2,1) model can be fit with:

arma_model <- arima(y, order = c(2, 0, 1))
summary(arma_model)

where the middle “0” indicates no differencing. For an ARIMA(1,1,1) model, which includes first differencing to remove trend, the code is:

arima_model <- arima(y, order = c(1, 1, 1))
summary(arima_model)

ccc

R packages

To conduct time series analysis we can use built in functions like ts() and decompose(). HoltWinters() also useful, now part of stats. Lots of specialized time series packages with advanced features, including forecast, timeSeries (Financial time series), season (Seasonal analysis of health data), and many others. For a more extensive package, see asta, which was designed to accompany the excellent book, now in its 5th edition, Time Series Analysis and its Applications: With R Examples, by Shumway and Stoffer (2025).

Time series data sets included in R and Rcmdr

R Code

data(co2, package="datasets")
co2 <- as.data.frame(co2)

#convert to time series data type with ts()
tCO2 <- ts(co2,frequency=12,start=c(1959),end=c(1997))
plot.ts(tCO2)

Figure 1. Time series plot of CO2 data set, the Keeling curve, from package datasets, comes with Rcmdr installation.

Other datasets included with R

carData::Arrests

carData::Bfox

carData::CanPop

Example

Get up-to-date CO2 data from NOAA as text file. Download to your computer, load and clean in your favorite spreadsheet app. Months came as numbers 1,2,3, etc., I changed to text, Jan, Feb, Mar, etc. I grabbed three columns: year, month, ppm for import to R.

head(maunaLoa)

R output

> head(maunaLoa)
  year month  ppm
1 1958 Mar 315.70
2 1958 Apr 317.45
3 1958 May 317.51
4 1958 Jun 317.24
5 1958 Jul 315.86
6 1958 Aug 314.93

However, it turns out the time series functions are easiest to work if only the ppm data are included.

tCO2 <- ts(maunaLoa[,"ppm"],frequency=12,start=c(1958,3),end=c(2020,10))
head(tCO2)

R output

> head(tCO2)
        Mar    Apr    May    Jun    Jul    Aug
1958 315.70 317.45 317.51 317.24 315.86 314.93

Get our plot (Figure 2).

plot(tCO2)

Figure 2. CO2 ppm monthly average data from NOAA, last data October 2020.

Seasonal time series comes with a trend component, a seasonal component, and a random component.

R code

dectCO2 <- decompose(tCO2)
head(dectCO2)
plot(dectCO2)

Figure 3. Observed (panel, top), trends over time (panel, second from top), seasonal changes (panel, second from bottom), and random error (panel, bottom).

Forecasting

Excellent resource at https://otexts.com/fpp2/

Exponential smoothing, weighted averages of past observations, weighted so that more recent observations are more influential.

Holt-Winters method extracts seasonal component (additive or multiplicative).

#set start value to value of first observation
tCO2cast <- HoltWinters(tCO2, l.start=315.42)

#Predict for next ten years. Because frequency in ts() was monthly, ten years is h=120
forecastCO2 <- forecast(tCO2cast, h=120)
plot(forecastCO2, fcol="red")

Figure 4. Data in black, predicted values in red (additive) shaded by confidence interval.

Questions

1. Write up three learning outcomes for this page. Hint: Point your favorite generative AI to this page and ask for help.
2. For the co2 dataset included in Rcmdr (co2, datasets), obtain forecast for year 2020 and compare against actual 2020 data (see Figure 2).
3. Positive clinical samples between September 2015 and November 2020 for flu virus in the USA are provided in the data set below (scroll or click here). The frequency of observations was weekly. Apply decompose() and obtain the seasonal and trend components of the data set. Which month does the peak positive sample occur?
4. Total pounds of fish (variable = Pounds) and pounds of Akule and Opelu (variable = Akule.Opelu) caught by commercial industry in Hawaiʻi, from 2000 to 2018 are provided in the data set below (scroll or click here). Apply decompose() and obtain the seasonal and trend components of the data set for Total pounds and again for Akule (Selar crumenophthalmus) and Opelu (Decapterus macarellus). Is there evidence for  trends, and if so, describe the trend. Is there evidence of seasonality? If so, which month did peak fishing occur?

Quiz Chapter 20.5

Time series

Start Quiz

Question

Confirm

✕

✓

Your answer:

Correct answer:

Next

Please wait..

↺

You got {{SCORE_CORRECT}} out of {{SCORE_TOTAL}}

Your Answers

References and suggested reading

Hamilton, J. D. (1994). Time Series Analysis. Princeton University Press.

Krebs, C. J., Boutin, S., Boonstra, R., Sinclair, A. R. E., Smith, J. N. M., Dale, M. R. T., Martin, K., & Turkington, R. (1995). Impact of food and predation on the snowshoe hare cycle. Science, 269(25 August), 1112–1115.

Kumar, A., Komaragiri, R., & Kumar, M. (2018). From Pacemaker to Wearable: Techniques for ECG Detection Systems. Journal of Medical Systems, 42(2), 34.

Chapter 6.4. Introduction to Time Series Analysis, in NIST/SEMATECH e-Handbook of Statistical Methods, https://www.itl.nist.gov/div898/handbook/index.htm, (first reviewed by Mike in 2019).

Mudelsee, M. (2016). Climate time series Classical statistics and bootstrap methods (2nd ed., Vol. 51). Springer.

Shumway, R. H., & Stoffer, D. S. (2025). Time Series Analysis and Its Applications: With R Examples (5th ed.). Springer.

Data set this page

Flu, extracted 28 Nov 2020 from https://gis.cdc.gov/grasp/fluview/fluportaldashboard.html

Data set in this page

Fish, Hawaiʻi state DLNR, Pounds refers to total catch, Akule.Opelu refers to pounds for the two kinds of fish.

/MD

Chapter 20 contents

• Additional topics
• Area under the curve
• Peak detection
• Baseline correction
• Surveys
• Time series
• Dimensional analysis
• Estimating population size
• Diversity indexes
• Survival analysis
• Growth equations and dose response calculations
• Plot a Newick tree
• Phylogenetically independent contrasts
• How to get the distances from a distance tree
• Binary classification
• Meta-analysis

/MD