如何在R语言中构建含线性成分的余弦模型（Cosinor Model）

Hi Katy! Welcome to Stack Overflow—no need to stress about etiquette gaps at all, we’re all here to learn and help each other out. Let’s walk through exactly how to build and fit your desired model step by step.

First, let's restate your model clearly for reference:

𝑦=𝐴cos{2𝜋((𝑡−∅)/τ)}+𝐶𝑡+𝐷
Fixed parameters: τ = 24 hours
Parameters to estimate: Amplitude (A), phase shift (∅), linear coefficient (C), constant term (D)

Step 1: Prepare Your Data #

First, make sure your time variable t is a numeric value representing hours (not just a datetime string). If you're working with datetime objects, you can convert them to hours since your starting point using the lubridate package:

# Install and load lubridate if you haven't already
# install.packages("lubridate")
library(lubridate)

# Example: Convert datetime to hours since the first observation
df$t_hours <- as.numeric(df$datetime - min(df$datetime), units = "hours")

(Replace df and datetime with your actual data frame and column names.)

Step 2: Define the Model Formula #

Since τ is fixed at 24, we can write the model directly in R syntax:

model_formula <- y ~ A * cos(2 * pi * (t_hours - phi)/24) + C * t_hours + D

Step 3: Set Initial Parameter Guesses #

Nonlinear models (like this one) need reasonable starting values to converge properly. Here’s how to make smart guesses:

• A (Amplitude): Use half the range of your y values: (max(df$y) - min(df$y))/2
• phi (Phase Shift): Guess the time of day where your y peaks (e.g., if peaks at 6 AM, use phi = 6)
• C (Linear Coefficient): Get the slope from a simple linear model: coef(lm(y ~ t_hours, data = df))[2]
• D (Constant Term): Calculate using the mean of y minus C * mean(t_hours)

Example code for initial values:

init_vals <- list(
  A = (max(df$y) - min(df$y))/2,
  phi = 6, # Adjust this based on your data's peak time
  C = coef(lm(y ~ t_hours, data = df))[2],
  D = mean(df$y) - coef(lm(y ~ t_hours, data = df))[2] * mean(df$t_hours)
)

Step 4: Fit the Model #

We’ll use R’s built-in nls() function for nonlinear least squares. If you run into convergence errors, the nlsLM() function from the minpack.lm package is more robust:

# Option 1: Basic nonlinear least squares
fit_nls <- nls(model_formula, data = df, start = init_vals)

# Option 2: Robust version (if nls fails)
# install.packages("minpack.lm")
library(minpack.lm)
fit_nlslm <- nlsLM(model_formula, data = df, start = init_vals)

Step 5: Inspect Results & Visualize #

To see your estimated parameters and model stats:

summary(fit_nls) # Or summary(fit_nlslm) if you used the robust version

Plot the fitted model against your raw data to check the fit:

# Install and load ggplot2 if needed
# install.packages("ggplot2")
library(ggplot2)

# Add predicted values to your data frame
df$predicted <- predict(fit_nls, newdata = df)

# Create the plot
ggplot(df, aes(x = t_hours, y = y)) +
  geom_point(alpha = 0.5, color = "gray") +
  geom_line(aes(y = predicted), color = "#2c3e50", linewidth = 1) +
  labs(x = "Time (Hours)", y = "Response Variable (y)", 
       title = "Fitted Cosine Model with Linear Trend") +
  theme_minimal()

Quick Tips #

• If the model fails to converge, tweak your initial guess for phi or A—small changes can make a big difference
• Compare this model to a simple linear model using anova(fit_nls, lm(y ~ t_hours, data = df)) to see if the cosine component adds significant value
• Ensure your time variable is continuous (avoid large gaps that might disrupt the periodicity assumption)

内容的提问来源于stack exchange，提问作者Katy J