Non-Normal Residuals in Real World Data I have a dataset that includes real world data (not experimental or survey data) for a set of countries year by year for 40 years. The data was collected by entities such as the World Bank and United Nations. We want to see how well this data predicts things like GDP. In the OLS case, R2=0.47 but the residuals are not anything close to normal. And bunch up on one side of the plot and have a strong linear trend as well--not random.
Curve fitting showed that a cubic distribution worked best. (A little effort yielded a lot of result.) R2 on the cubic model was .59. SPSS only gives unstandardized residuals for nonlinear regression (Maybe because standardized doesn't make sense?) Plotting the unstandardized residuals is still not normal. Shapiro-Wilk is still <.001. The Q-Q plot kinda follows the expected line though, see photo. What should I do about this?


One option since you are considering nonlinear regression is a generalized additive model (GAM). I have provided an example below. First I loaded the mgcv package for running GAMs and MASS for the mcycle dataset, which is referenced in my link at the end of this answer.
#### Libraries ####
library(mgcv)
library(MASS)

I then fit a model. Its fairly similar to a normal regression fit using lm with some minor tweaks. First, I added a simple spline term here (basically a function which creates a smooth function for nonlinear modeling), then set the number of knots to 10, basically stating that the maximum amount of curving points should be around this number. Thereafter I just add the dataset and set the method to either REML or ML (these are the two suggested methods to use).
#### Fit Model ####
fit <- gam(
  accel ~ s(times,
            k=10),
  data = mcycle,
  method = "REML"
)

From there you can just summarize the model like an lm object.
#### Summarize and Check ####
summary(fit)

Which gives you these terms:
Family: gaussian 
Link function: identity 

Formula:
accel ~ s(times, k = 10)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -25.546      1.951  -13.09   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
           edf Ref.df    F p-value    
s(times) 8.625  8.958 53.4  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.783   Deviance explained = 79.7%
-REML = 616.14  Scale est. = 506.35    n = 133

You'll see the summary looks different. Here, if you chose linear terms, they would be listed with the intercept under the parametric section. We only have one slope term, the smooth fit function for times, so we skip to that in the Approximate Smoothness section. Here we see the term is significant and the edf is around 8.625, indicating a fair amount of curvature. We can visualize it by simply plotting it directly. Here I've added the residual argument to show the data points and changed the point shapes.
#### Plot Model ####
plot(fit,
     residuals = T,
     pch=21)


You can see that the data clearly has a nonlinear function, so fitting a linear model would be fairly bad. We can also check to make sure our model doesn't need to be adjusted (this is simply a heuristic check and if p values are flagged you can consider adjusting your splines to fit better). It automatically spits out a plot of the fitted/response values as well.
gam.check(fit)

The summary and plot are seen below:

Method: REML   Optimizer: outer newton
full convergence after 7 iterations.
Gradient range [-1.436342e-06,1.179834e-06]
(score 616.142 & scale 506.3529).
Hessian positive definite, eigenvalue range [3.329028,65.73378].
Model rank =  10 / 10 

Basis dimension (k) checking results. Low p-value (k-index<1) may
indicate that k is too low, especially if edf is close to k'.

           k'  edf k-index p-value
s(times) 9.00 8.62    1.15    0.95

You should get 4 plots, and the qq plot this time looks normal.

There appear to be no issues with our check, so we are free to move on.
Useful Resources
For learning all about GAMs, the DataCamp course is great, but this is also an excellent free resource:
https://noamross.github.io/gams-in-r-course/
