$$\begin{tabular}{|1|1|}\hline Year & Population\\ \hline 1959&4835\\ 1960&4970\\ 1961&5085\\ 1962&5160\\ 1963&5310\\ 1964&5260\\ 1965&5235\\ 1966&5255\\ 1967&5235\\ 1968&5210\\ 1969&5175\\ \hline \end{tabular}$$
First we import the data into R:
Year <- c(1959, 1960, 1961, 1962, 1963, 1964, 1965, 1966, 1967, 1968, 1969)
Population <- c(4835, 4970, 5085, 5160, 5310, 5260, 5235, 5255, 5235, 5210, 5175)
Now we create the dataframe named sample1:
sample1 <- data.frame(Year, Population)
sample1
Year Population
1 1959 4835
2 1960 4970
3 1961 5085
4 1962 5160
5 1963 5310
6 1964 5260
7 1965 5235
8 1966 5255
9 1967 5235
10 1968 5210
11 1969 5175
At this point may be useful to chart these values, to observe the trend and take an idea of the final polynomial function. For convenience we modify the column
Year, creating a neighborhood of zero, thus:
sample1$Year <- sample1$Year - 1964
sample1
Year Population
1 -5 4835
2 -4 4970
3 -3 5085
4 -2 5160
5 -1 5310
6 0 5260
7 1 5235
8 2 5255
9 3 5235
10 4 5210
11 5 5175
Put the values on a chart
plot(sample1$Year, sample1$Population, type="b")
At this point we can start with the search for a polynomial model that adequately approximates our data. First, we specify that we want a polynomial function of X, ie a raw polynomial , is different from the orthogonal polynomial. This is an important addition because the controls and the results will change in the two cases R. So we want a function of X like:
$$f(x)=\beta_0+\beta_1x+\beta_2x^2+\beta_3x^3+ ... +\beta_nx^n$$
At what degree of the polynomial stop? Depends on the degree of precision that we seek. The greater the degree of the polynomial, the greater the accuracy of the model, but the greater the difficulty in calculating; we must also verify the significance of coefficients that are found. But let's get straight to the point.
In R for fitting a polynomial regression model (not orthogonal), there are two methods, among them identical. Suppose we seek the values of beta coefficients for a polynomial of degree 1, then 2nd degree, and 3rd degree:
fit1 <- lm(sample1$Population ~ sample1$Year)
fit2 <- lm(sample1$Population ~ sample1$Year + I(sample1$Year^2))
fit3 <- lm(sample1$Population ~ sample1$Year + I(sample1$Year^2) + I(sample1$Year^3))
Or we can write more quickly, for polynomials of degree 2 and 3:
fit2b <- lm(sample1$Population ~ poly(sample1$Year, 2, raw=TRUE))
fit3b <- lm(sample1$Population ~ poly(sample1$Year, 3, raw=TRUE))
The function
poly is useful if you want to get a polynomial of high degree, because it avoids explicitly write the formula. If we specify raw=TRUE, the two methods provide the same output, but if we do not specify raw=TRUE (or rgb(153, 0, 0);">raw=F), the function poly give us the values of the beta parameters of an orthogonal polynomials, which is different from the general formula I wrote above, although the models are both effective.Let's look at the output.
summary(fit2)
## or summary(fit2b)
Call:
lm(formula = sample1$Population ~ sample1$Year + I(sample1$Year^2))
Residuals:
Min 1Q Median 3Q Max
-46.888 -18.834 -3.159 2.040 86.748
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5263.159 17.655 298.110 < 2e-16 ***
sample1$Year 29.318 3.696 7.933 4.64e-05 ***
I(sample1$Year^2) -10.589 1.323 -8.002 4.36e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 38.76 on 8 degrees of freedom
Multiple R-squared: 0.9407, Adjusted R-squared: 0.9259
F-statistic: 63.48 on 2 and 8 DF, p-value: 1.235e-05
The output of
summary(fit2b) is the same. We obtained the values of beta0 (5263,159), beta1 (29,318) and beta2 (-10,589), which appear to be significant AII 3. The equation of polynomial of degree 2 of our model is:$$f(x)=5263.1597+29.318x-10.589x^2$$
If we want a polynomial of 3rd degree, we have:
summary(fit3)
## of summary(fit3b)
Call:
lm(formula = sample1$Population ~ sample1$Year + I(sample1$Year^2) +
I(sample1$Year^3))
Residuals:
Min 1Q Median 3Q Max
-32.774 -14.802 -1.253 3.199 72.634
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5263.1585 15.0667 349.324 4.16e-16 ***
sample1$Year 14.3638 8.1282 1.767 0.1205
I(sample1$Year^2) -10.5886 1.1293 -9.376 3.27e-05 ***
I(sample1$Year^3) 0.8401 0.4209 1.996 0.0861 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 33.08 on 7 degrees of freedom
Multiple R-squared: 0.9622, Adjusted R-squared: 0.946
F-statistic: 59.44 on 3 and 7 DF, p-value: 2.403e-05
The equation is:
$$f(x)=5263.1585+14.3638x-10.5886x^2+0.8401x^3$$
In the latter case, however, the coefficients beta1 and beta3 are not significant, then the best model is the polynomial of 2nd degree. Furthermore look at the Multiple R-squared: in the 2nd degree model it is 94.07%, while in the 3rd degree model it is 96.22%. It seems that there has been an increase in accuracy of the model, but it is a significant increase? We can compare the two model with an ANOVA table:
anova(fit2, fit3)
Analysis of Variance Table
Model 1: sample1$Population ~ sample1$Year + I(sample1$Year^2)
Model 2: sample1$Population ~ sample1$Year + I(sample1$Year^2) + I(sample1$Year^3)
Res.Df RSS Df Sum of Sq F Pr(>F)
1 8 12019.8
2 7 7659.5 1 4360.3 3.9848 0.0861 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Since the p-value is greater than 0.05, we accept the null hypothesis: there wasn't a significant improvement of the model.
The biggest problem now is to represent graphically the result. In fact, R does not exist (as far as I know) a function for plotting polynomials found. We must therefore proceed with graphic artifacts still valid, but somewhat laborious.
First, we plotted the values, with the command seen before. This time only display the lines and not points, for convenience graphics:
plot(sample1$Year, sample1$Population, type="l", lwd=3)
Now add to this chart the progress of the 2nd degree polynomial, in this way:
points(sample1$Year, predict(fit2), type="l", col="red", lwd=2)
The function
predict() compute the Y values given the X values. The the coordinates are linked with lines. Is not plotted the continuous, but the discrete. With a few values, this method is highly debilitating.
Let's add the graph of the polynomial of 3rd degree:
points(sample1&Year, predict(fit3), type="l", col="blue", lwd=2)
As you can see the two models have very similar trends.
If we would instead obtain the graph of continuous functions obtained, we proceed in this manner. First you create the polynomial equation we previously found:
pol2 <- function(x) fit2$coefficient[3]*x^2 + fit2$coefficient[2]*x + fit2$coefficient[1]
Remember that:
- coefficient[1] = beta0
- coefficient[2] = beta1
- coefficient[3] = beta2
and so on.
At this point we plotted the coordinates of sample1 and then the created curve with
curve(x):
plot(sample1$Year, sample1$Population, type="p", lwd=3)
pol2 <- function(x) fit2$coefficient[3]*x^2 + fit2$coefficient[2]*x + fit2$coefficient[1]
curve(pol2, col="red", lwd=2)
The point, however, disappear, but we can replace them with the command
points:
points(sample1$Year, sample1$Population, type="p", lwd=3)
A note: you must follow the order of commands as I have described, otherwise the function
curve creates a wrong graph. So summarizing the commands to get the continuous function, and the experimental points on the same graph are the following:
plot(sample1$Year, sample1$Population, type="p", lwd=3)
pol2 <- function(x) fit2$coefficient[3]*x^2 + fit2$coefficient[2]*x + fit2$coefficient[1]
curve(pol2, col="red", lwd=2)
points(sample1$Year, sample1$Population, type="p", lwd=3)
The graph we get is the following:
Now draw the graph of the polynomial of 3rd degree:
plot(sample1$Year, sample1$Population, type="p", lwd=3)
pol3 <- function(x) fit3$coefficient[4]*x^3 + fit3$coefficient[3]*x^2 + fit3$coefficient[2]*x + fit3$coefficient[1]
curve(pol3, col="red", lwd=2)
points(sample1$Year, sample1$Population, type="p", lwd=3)
