# Weighted least squares, checking normality

• Assuming I want to perform a linear regression, I can do:
# Generate some example data
set.seed(123)
x <- rnorm(100, mean = 5, sd = 2)
y <- 3 + 2 * x + rnorm(100, mean = 0, sd = 1)

# Fit the linear regression model
model <- lm(y ~ x)

# Plot the data with the regression line
plot(x, y, main = "Scatter Plot with Regression Line", xlab = "X", ylab = "Y")
abline(model, col = "red")


# Obtain the residuals
residuals <- model$residuals # Check the normality of residuals with a Q-Q plot qqnorm(residuals) qqline(residuals, col = "red") # Check the normality of residuals with the Shapiro-Wilk test shapiro.test(residuals)   Shapiro-Wilk normality test data: residuals W = 0.9748, p-value = 0.05204  Here we assume the data have a common variance (unknown). In the second case, I will assume the error for each predictor is estimated. • Now let's assume, I have an experiment where I collect for each x value (e.g. week 1 week 2, week 3), some data (y). set.seed(0) # Generate random normal values y1 <- rnorm(10, mean = 5, sd = 2) y2 <- rnorm(12, mean = 8, sd = 2) y3 <- rnorm(14, mean = 12, sd = 2) y_values <- c(y1, y2, y3) x_values <- c(rep(1, length(y1)), rep(2, length(y2)), rep(3, length(y3))) # Plot the values plot(x_values, y_values, xlab="x", ylab="y", pch=19, col=as.factor(x_values))  I want to see if there is a linear relation between the means of y_i and the x_i. As I have multiple groups of data (e.g., y1, y2, y3) and I want to check the normality before performing a linear regression, should I treat each group separately ? That is : # Check normality for each group shapiro.test(y1) shapiro.test(y2) shapiro.test(y3)  > shapiro.test(y1) Shapiro-Wilk normality test data: y1 W = 0.96752, p-value = 0.867 > shapiro.test(y2) Shapiro-Wilk normality test data: y2 W = 0.95233, p-value = 0.6713 > shapiro.test(y3) Shapiro-Wilk normality test data: y3 W = 0.96314, p-value = 0.7743  # Plot Q-Q plots for each group par(mfrow = c(1, 3)) qqnorm(y1); qqline(y1, col = "red") qqnorm(y2); qqline(y2, col = "red") qqnorm(y3); qqline(y3, col = "red") par(mfrow = c(1, 1))  And then, how do I check for normality when performing ? # Calculate the means and SEM for each group mean_y1 <- mean(y1); sem_y1 <- sd(y1) / sqrt(length(y1)) mean_y2 <- mean(y2); sem_y2 <- sd(y2) / sqrt(length(y2)) mean_y3 <- mean(y3); sem_y3 <- sd(y3) / sqrt(length(y3)) # Create a data frame for plotting and modeling data <- data.frame( x = 1:3, mean = c(mean_y1, mean_y2, mean_y3), sem = c(sem_y1, sem_y2, sem_y3) ) # Fit the nonlinear model using nls nlmodel <- nls(mean ~ A * x + B, data = data, start = list(A = 1, B = 1), weights = 1/sem^2) summary(nlmodel) # Create a function to predict values using the fitted model predict_nlmodel <- function(A, B, newdata) { A * newdata + B } # Get the fitted parameters params <- coef(nlmodel) # Plot the data and the fitted nonlinear model ggplot(data, aes(x = x, y = mean)) + geom_point(size = 3) + geom_errorbar(aes(ymin = mean - sem, ymax = mean + sem), width = 0.1) + stat_function(fun = predict_nlmodel, args = list(A = params["A"], B = params["B"]), color = "blue") + labs(title = "Nonlinear Weighted Regression", x = "X", y = "Mean Value") + theme_minimal()  (I use nls but the model is a linear model here). Should I do something like: # Extract residuals and fitted values residuals <- residuals(nlmodel) fitted <- fitted(nlmodel) # Q-Q Plot ggplot(data.frame(sample = residuals), aes(sample = sample)) + stat_qq() + stat_qq_line() + labs(title = "Q-Q Plot of Residuals", x = "Theoretical Quantiles", y = "Sample Quantiles") + theme_minimal()  Which I guess is wrong cause each predictor does not have the same variance. • The iid. normal assumption applies to your error term, not the data, so it's only the residuas you should be concerned about. Even then formally testing for normality is not recommended; your first example captures the data-generating mechanism exactly and still gives you$P=0.052$. Commented Jul 27 at 14:46 • Thanks. So the second approach is wrong ? Commented Jul 27 at 14:48 • It doesn't make much sense, no. A standard linear model can capture shifts in mean as in your second example (not clear if these 'groups' are categories or continuous but the way you've set it up either works). The assumption remains that these groups have the same conditional variance, also still met here. If that weren't the case it would start showing as systematic trends in your residuals, which is evidence your model is misspecified. Commented Jul 27 at 15:01 • In experimental physics or biology, I guess that this what people do. They collect data (y) for different x values and they may want to see if there is a linear relationship between the means of y_i and x_i. To do so they calculate the mean for each y_i, calculate the standard error on the mean and fit. Then, they calculate a reduced chi2 to access the goodness of the fit. If the reduced chi2 is close to 1, this means the linear model is capable of capturing experimental trends. But in this case, I am not sure how one should check for the normality. Commented Jul 27 at 15:11 • Welcome to CV! The correct residuals for nlmodel are determined by the original data, not the residuals based on the means. Your question amounts to testing the hypothesis that each group of residuals has a common variance versus the alternative that at least one group has a different variance. Normality testing can be done but is rarely useful or even correctly interpreted--you can find many threads here on CV about that topic. – whuber Commented Jul 27 at 16:26 ## 1 Answer Let's take, for simplicity, the case of a simple OLs regression. The model is $$y_i=b_0+b_1.x_i+\epsilon_i$$, or $$\hat y_i=b_0+b_1.x_i$$. You can perform a regression, and compute the coefficients $$b_i$$, and compute the $$r$$ or $$R^2$$ etc., all you want, without any distributional assumption. All these coefficients are descriptive statistics, such as the mean $$\bar x$$, standard deviation $$s$$, median, skewness, etc. They can be validly computed w/o assuming normality of any distribution. All these coefficients ($$b_i, r$$, ...) are just the result of applying a mathematical formula to the data, just like for $$\bar x$$,$$s$$, etc. And like the mean is the point which minimizes the sum of squared errors, the regression line minimizes the sum of squared errors between the predicted $$\hat y_i$$'s and the observed $$y_i$$'s. Again, no normality of anything. If you are concerned by higher order linear models, just replace point (mean, for 1D) and line (regression line, for 2D), by plane (3D), or hyperplane (nD) in the above paragraph. Now, if you want to make inferences about these coefficients, e.g. Confidence Intervals (CI's) around the coefficients, or p-values (CI's and p's are just flip sides of the same coin), now we do need to make distributional assumptions (sort of...). In this case, it is the residuals ($$\hat y_i-y_i$$) which are assumed to be normaly distributed. But not the $$y_i's$$, as you did. Note that you can (and often do) run a regression where there is only 1 value of $$y_i$$ at each $$x_i$$; how would we test normality of the y's at each x, as you did? So, only normality of the residuals, i.e. $$(\hat y_i-y_i)'s$$. And even if the residuals are not very normal looking, the CI's, and p-values will not be so wrong as to horribly mislead you. In your example with the 3 $$x_i's$$, from the graph alone, it should be clear that the slope is not the result or random noise, and that factor is significant. You may want to look at this excellent response on CV, which goes over the various regression assumptions. Bottom line; OLS will still give you a very practical usable answer. • Thanks. But I still do not understand what I should do for the second model ... In the first model, when you have 1 value of 𝑦𝑖 at each 𝑥, I guess people are checking the normality on all yi as I did it. As for the second model, I thought all I should check the normality at each x_i value (stats.stackexchange.com/questions/148803/…). Commented Jul 27 at 17:03 • You would do exactly the same thing as fort the 1st case; you combine all the residuals together; compute all$(\hat y_i - y_i)$'s, and check normality. Note that$x_i$does not appear here; so no difference if each$y_i$has a unique$x_i$, or if one or more$y_i$share the same$x\$... Commented Jul 28 at 1:30
• Ok. Thanks. I think I understood: Ordinary Least Squares (OLS) regression does not require the dependent variable or the residuals to be normally distributed for estimating the coefficients. However, the assumption of normality becomes important when you want to determine Confidence Intervals. (Weighted) WLS also does not require the dependent variable or the residuals to be normally distributed for estimating the coefficients. OLS estimation is equivalent to MLE when the errors are normally distributed. When there is heteroscedasticity in the residuals, one should perform a WLS. Commented Jul 28 at 11:11
• But I still do not understand you have to combine all residuals together in case 2. Looking at stats.stackexchange.com/questions/143705/…, I would say that the normality assumption is for every y_i. I guess (but I may be wrong) that is kind if 'not correct' to aggregate when having heteroscedasticity. Commented Jul 28 at 11:18
• @Laut567, you check normality by of all residuals combined, regardless of homoscedasticity or not. You check all the residuals, at each predictor variable for homoscedasticity, regardless of normality or not. But before you spend too much time validating these assumptions, you should read the following 2 blogs; they will give you some good insights about how important the assumptions are (or are not), and what happens when they are not met; janhove.github.io/posts/2019-04-11-assumptions-relevance and economictheoryblog.com/2015/04/01/ols_assumptions. Commented Jul 31 at 19:43