线性回归建模物种迁徙：如何同时满足假设并修正自相关？

Hey there! Let's work through this problem step by step—you're dealing with proportional citizen science data, linear regression assumptions, and autocorrelation from monthly time series, so we need to address each of these to build a reliable model.

Your response variable is a percentage (derived from monthly sightings), which violates the linear regression assumption of normally distributed residuals (since percentages are bounded between 0-100). Here's what to do:

• Convert to 0-1 proportions: First, scale your percentage down to a 0-1 range (divide by 100).
• Choose the right model type:
  • If your proportions have no exact 0s or 1s, use beta regression (via the betareg package). Beta distributions are designed specifically for proportional data, so this will align better with your response variable's distribution than a standard lm().
  • If you do have 0s or 1s, adjust the values slightly (e.g., add 0.5 / total_monthly_sightings to 0s, subtract the same from 1s) before using beta regression, or use a beta-binomial regression to account for extra variability from citizen science data.

Since your data is monthly, autocorrelation (residuals being correlated across time) is almost guaranteed, which breaks the lm() assumption of independent errors. Here are your best fixes:

Option 1: Generalized Least Squares (GLS) with Autocorrelation Structure #

Use the nlme package's gls() function to explicitly model the time-dependent correlation. For example, if you suspect first-order autocorrelation (AR(1), where each month's value correlates with the previous one):

library(nlme)
# First, logit-transform your proportion if using GLS (to approximate normality)
your_data$logit_prop <- log(your_data$proportion / (1 - your_data$proportion))

# Fit GLS with AR(1) correlation
model_gls <- gls(logit_prop ~ SST + Chl_a, 
                 data = your_data,
                 correlation = corAR1(form = ~ month)) # 'month' is your time variable

This model will estimate the autocorrelation parameter and adjust the model to account for it, making your residuals independent.

Option 2: Newey-West Standard Errors (Quick Fix for Linear Models) #

If you want to stick with a linear model (e.g., after logit transformation), use Newey-West standard errors to correct for autocorrelation and heteroscedasticity without changing the model coefficients. Use the sandwich package:

library(sandwich)
library(lmtest)

model_lm <- lm(logit_prop ~ SST + Chl_a, data = your_data)
# Calculate Newey-West SEs
nw_se <- neweywest(model_lm, lag = 1) # Adjust lag based on your ACF plot
# Get corrected p-values
coeftest(model_lm, vcov = nw_se)

This is a good option if you just need reliable statistical inference, not necessarily to fix the residual independence assumption for prediction.

Option 3: ARIMAX Models (For Time Series Forecasting + Covariates) #

If your goal includes predicting future migration patterns, use an ARIMAX model (ARIMA with exogenous covariates like SST and Chl-a). The forecast package makes this easy:

library(forecast)
# Logit-transform the proportion first
model_arimax <- auto.arima(your_data$logit_prop, 
                           xreg = your_data[, c("SST", "Chl_a")],
                           seasonal = TRUE) # Add if you see seasonal patterns

This model automatically selects the best ARIMA structure to handle autocorrelation while including your environmental predictors.

Once you've chosen a model, you need to confirm assumptions are met:

• For GLS/beta regression:
  • Linearity: Plot transformed response vs. each covariate (e.g., plot(logit_prop ~ SST, data = your_data)) to check for linear relationships. If you see curves, add polynomial terms (e.g., poly(SST, 2)).
  • Residual independence: Use acf(resid(model_gls)) to check that residuals have no remaining autocorrelation.
  • Homoscedasticity: Plot residuals vs. fitted values—no funnel shapes should be present. Use the Breusch-Pagan test if needed.
  • Distribution: For beta regression, check that residuals follow a beta distribution (use Q-Q plots tailored to beta data, or package-specific diagnostic functions).

1. Preprocess your proportion data (scale to 0-1, adjust 0s/1s if needed).
2. Explore time trends and covariate relationships with plots.
3. Choose a model that addresses both proportional data and autocorrelation (GLS with AR(1) or beta regression + Newey-West, or ARIMAX).
4. Diagnose residuals to ensure all assumptions are satisfied.
5. Interpret coefficients: For logit-transformed models, exponentiate coefficients to get odds ratios; for beta regression, interpret via the link function you used.

内容的提问来源于stack exchange，提问作者Jo Harris