Great question—this is a common confusion because Cohen's d is typically taught for t-tests or ANOVA scenarios, but extending it to multiple (including polynomial) regression requires clarifying exactly which effect you want to quantify. The effsize package doesn't work with your regression formula because it's designed for group comparisons, not continuous predictor models with polynomial terms. Here are two robust approaches tailored to your model:

1. Use Standardized Regression Coefficients (Beta) as a Cohen's d Analog #

For continuous predictors, a standardized regression coefficient (beta) directly tells you how many standard deviations the dependent variable (mwb) changes when the predictor (watch_wd) increases by one standard deviation—this aligns closely with the intuition of Cohen's d for continuous effects.

To calculate this:

# Standardize your dependent and predictor variables
new_data_std <- as.data.frame(scale(new_data[, c("mwb", "watch_wd")]))

# Fit the polynomial regression on standardized data
lr_watchweek_std <- lm(mwb ~ watch_wd + I(watch_wd^2), data = new_data_std)

# View standardized coefficients (beta values)
summary(lr_watchweek_std)$coefficients

• The coefficient for watch_wd is the linear effect's standardized size (e.g., a value of 0.3 means mwb increases by 0.3 SDs per 1 SD increase in watch_wd).
• The coefficient for I(watch_wd^2) quantifies the standardized quadratic effect.

2. Calculate Cohen's d for Key Predictor Values #

If you want a single Cohen's d that summarizes the overall effect of watch_wd (including its quadratic term), you can compare predicted mwb values at meaningful extremes of watch_wd (e.g., ±1 standard deviation from the mean) and divide that difference by the SD of mwb:

# Calculate mean and SD for watch_wd and mwb
wd_mean <- mean(new_data$watch_wd, na.rm = TRUE)
wd_sd <- sd(new_data$watch_wd, na.rm = TRUE)
mwb_sd <- sd(new_data$mwb, na.rm = TRUE)

# Define extreme values of watch_wd (±1 SD)
low_wd <- wd_mean - wd_sd
high_wd <- wd_mean + wd_sd

# Predict mwb at these extremes
pred_low <- predict(lr_watchweek, newdata = data.frame(watch_wd = low_wd))
pred_high <- predict(lr_watchweek, newdata = data.frame(watch_wd = high_wd))

# Compute Cohen's d
cohen_d <- (pred_high - pred_low) / mwb_sd
cat("Cohen's d for watch_wd (±1 SD):", cohen_d, "\n")

This d value represents how many standard deviations mwb changes when watch_wd moves from one standard deviation below the mean to one above, accounting for both linear and quadratic effects.

Why effsize Didn't Work #

The cohen.d() function in effsize is built for group-based comparisons (e.g., comparing mwb between two watch_wd groups). Your model uses a continuous predictor with a polynomial term, so there's no inherent "group" to compare—hence the incompatibility.

内容的提问来源于stack exchange，提问作者Sarah.d