17.4 – OLS, RMA, and smoothing functions

Introduction

OLS or ordinary least squares is the most commonly used estimation procedure for fitting a line to the data. For both simple and multiple regression, OLS works by minimizing the sum of the squared residuals. OLS is appropriate when the linear regression assumptions LINE apply. In addition, further restrictions apply to OLS including that the predictor variables are fixed and without error. OLS is appropriate when the goal of the analysis is to retrieve a predictive model. OLS describes an asymmetric association between the predictor and the response variable: the slope b_X for Y ~ b_XX will generally not be the same as the slope b_Y for X ~ b_YY.

OLS is appropriate for assessing functional relationships (i.e., inference about the coefficients) as long as the assumptions hold. In some literature, OLS is referred to as a Model I regression.

Generalized Least Squares

Generalized linear regression is an estimation procedure related to OLS but can be used either when variances are unequal or multicollinearity is present among the error terms.

Weighted Least Squares

A conceptually straightforward extension of OLS can be made to account for situation where the variances in the error terms are not equal. If the variance of Yi varies for each Xi, then a weighting function based on the reciprocal of the estimated variance may be used.

then, instead of minimizing the squared residuals as in OLS, the regression equation estimates in weighted least squares minimizes the squared residuals summed over the weights.

Weighted least squares is a form of generalized least squares. In order to estimate w_i, however, multiple values of Y for each observed X must be available.

Reduced Major Axis

There are many alternative methods available when OLS may not be justified. These approaches, collectively, may be called Model II regression methods. These methods are invariably invoked in situations in which both Y and X variables have random error associated with them. In other words, the OLS assumption that the predictor variables are measured without error is violated. Among the more common methods is one called Reduced Major Axis or RMA.

Smoothing functions

data set: atmospheric carbon dioxide (CO₂) readings Mauna Loa. Source: https://gml.noaa.gov/ccgg/trends/data.html

Fit curves without applying a known formula. This technique is called smoothing and, while there are several versions, the technique involves taking information from groups of observations — weighted averaged — and using these groups to estimate how the response variable changes with values of the independent variable. Smoothing is used to help reveal patterns, to emphasize trends by reducing noise — clearly, caution needs to be employed as smoothing necessarily hides outlier data, which can themselves be important. Smoothing techniques by name include kernel, loess, and spline. Default in the scatter plot command is loess.

CO₂ in parts per million (ppm) plotted by year from 1958 to 2014 the first CO₂ readings were recorded in April 1958 (Fig 1).

A few words of explanation for Figure 1. The green line shows the OLS line, the red line shows the loess smoothing with a smoothing parameter of 0.5 (in Rcmdr the slider value reads “5”).

Figure 1. CO₂ in parts per million (ppm) plotted by year from 1958 to 2014

R command was started with option settings available in Rcmdr context menu for scatterplot, then additional commands were added

scatterplot(CO2~Year, reg.line=lm, grid=FALSE, smooth=TRUE, spread=FALSE, boxplots=FALSE, span=0.05, lwd=2, xlab="Months since 1958", ylab="CO2 ppm", main="CO2 at Mauna Loa Observatory, April 1958 - April 2014", cex=1, cex.axis=1.2, cex.lab=1.2, pch=c(16), data=StkCO2MLO)

The next plot is for ppm CO₂ by month for the year 2013. The plot shows the annual cycle of atmospheric CO₂ in the northern hemisphere.

Figure 2. Plot of ppm CO₂ by month for the year 2013.

Again, the smoothing parameter was set to 0.5 and the loess function is plotted in red (Fig 2).

Loess is an acronym short for local regression. Loess is a weighted least squares approach which is used to fit linear or quadratic functions of the predictors at the centers of neighborhoods. The radius of each neighborhood is chosen so that the neighborhood contains a specified percentage of the data points. This percentage of data points is referred to as the smoothing parameter and this parameter may differ for different neighborhoods of points. The idea of loess, in fact, any smoothing algorithm, is to reveal pattern within a noisy sequence of observations. The smoothing parameter can be set to different values, between 0 and 1 is typical.

To get a sense of what the parameter does, Figure 3 takes the same data as in Figure 2, but with different values of the smoothing parameter (Fig 3).

Figure 3. Plot with different smoothing values (0.5 to 10.0).

plot key:

The R code used to generate Figure 3 plot was

spanList = c(0.5, 0.75, 1, 2, 10)
reg1 = lm(ppm~Month)
png(filename = "RplotCO2mo.png", width = 400, height = 400, units = "px", pointsize = 12, bg = "white")
plot(Month,ppm, cex=1.2, cex.axis=1.2, cex.lab=1.2, pch=c(16), xlab="Months", ylab="CO2 ppm", main="CO2 levels 2013")
abline(reg1,lwd=2,col="green") 
for (i in 1:length(spanList))
{
ppm.loess <- loess(ppm~Month, span=spanList[i], Dataset)  
ppm.predict <- predict(ppm.loess, Month)
lines(Month,ppm.predict,lwd=2,col=i)
}

The CO2 data constitutes a time series. Instead of loess, a simple moving average would be a more natural way to reveal trends. In principle, take a set of nearby points (odd number of points best, keeps the calculation symmetric) and calculate the average. Next, shift the points by a specified time interval (eg, 7 days), and recalculate the average for the new set of points. See Chapter 20.5 for Time series analysis.

Questions

1. Write up three learning outcomes for this page. Hint: Point your favorite generative AI to this page and ask for help.
2. This is bio-statistics, so I gotta ask: What ecological/environmental process explains the shape of the relationship between ppm CO₂ and months of the year as shown in Figure 2? Hint: NOAA Global Monitoring Laboratory responsible for the CO2 data is located at Mauna Loa Observatory, Hawaiʻi (lat: 19.52291, lon: -155.61586).
3. As I wrote this question (January 2022), we were 22 months since WHO declared Covid-19 a pandemic (CDC timeline). Omicron variant is now dominant; Daily case counts State of Hawaiʻi from 1 November 2021 to 15 January 2022 reported in data set table.
   1. Make plot like Figure 2 (days instead of months)
   2. Apply different loess smoothing parameters and re-plot the data. Observe and describe the change to the trend between case reports and days.

Quiz Chapter 17.4

OLS, RMA, and smoothing functions

Start Quiz

Question

Confirm

✕

✓

Your answer:

Correct answer:

Next

Please wait..

↺

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

Your Answers

Data set

Covid-19 cases reported State of Hawaiʻi from 1 November 2021 to 15 January 2022 (data extracted from Wikipedia)

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
• ANCOVA – analysis of covariance
• 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