This is an assignment about Statistics course by using the R Studio.

Chapter 16: Inference for Regression

Climate Change

The earth has been getting warmer. Most climate scientists agree that one important cause of the warming is the increase in atmospheric levels of carbon dioxide (CO2), a green house gas. Here is part of a regression analysis of the mean annual CO2 concentration (CO2) in the atmosphere, measured in parts per thousand (ppt), at the top of Mauna Loa in Hawaii and the mean annual air temperature (Temp) over both land and sea across the globe, in degrees Celsius.

Let’s first read the dataset into R climate <- read.table('Climate_Change.txt', sep = '\t', header = TRUE)

and take a look at the data structure: str(climate)

## 'data.frame': 29 obs. of 3 variables: ## $ year: int 1980 1981 1982 1983 1984 1985 1986 1987 1988 1989 ... ## $ Temp: num 14.2 14.3 14.1 14.3 14.1 ... ## $ CO2 : num 339 340 341 342 344 ...

We see three variables, which are year, Temp (mean annual air temperature) and CO2 (mean annual CO2 concentration), and there are 29 observations in each variable.

We now take Temp as the response variable and CO2 the predictor variable, to study their relationship. To see if linear regression is appropriate, we make a scatterplot of Temp against CO2 plot(climate$CO2, climate$Temp, xlab = 'CO2 Concentration', ylab = 'Temperature')

340 350 360 370 380

1 4

.1 1

4 .3

1 4

.5

CO2 Concentration

Te m

p e

ra tu

re

It seems reasonable to fit a linear model to the dataset, because both variables are quantitative, the data points show a linear pattern, and there is no outlier. So, let’s fit the model: imod <- lm(Temp ~ CO2, data = climate)

1

The summary of the fitted model is given by summary(imod)

## ## Call: ## lm(formula = Temp ~ CO2, data = climate) ## ## Residuals: ## Min 1Q Median 3Q Max ## -0.16809 -0.07972 0.00194 0.07013 0.18532 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 10.707076 0.481006 22.260 < 2e-16 *** ## CO2 0.010062 0.001336 7.534 4.19e-08 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 0.09847 on 27 degrees of freedom ## Multiple R-squared: 0.6776, Adjusted R-squared: 0.6657 ## F-statistic: 56.76 on 1 and 27 DF, p-value: 4.192e-08

which contains a lot of information. We see that R2 = 0.6776 and the SD of residuals se = 0.09847 (the estimator of population standard deviation σ) with 27 degrees of freedom. In Coefficients section we see the intercept b0 = 10.71 and the slope b1 = 0.01. Their standard errors are SE(b0) = 0.481 and SE(b1) = 0.00134. Their t-test statistics are t0 = b0/SE(b0) = 22.26 and t1 = b1/SE(b1) = 7.534. Their corresponding (two-tailed) p-values are very small (<2e-16 and 4.19e-08). As a result, we reject H0 : β1 = 0 and conclude there is a positive correlation between Temp and CO2. The b1 = 0.01 can be interpreted as follows: The air temperature will increase by 0.01 degrees Celsius on average if the CO2 concentration in the atmosphere increases by 1 ppt. If only focus on the coefficients we may do summary(imod)$coefficients

## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 10.70707624 0.48100637 22.259739 6.641063e-19 ## CO2 0.01006241 0.00133563 7.533828 4.191615e-08

To obtain the confidence interval for each regression coefficient, we may do confint(imod, level = 0.95)

## 2.5 % 97.5 % ## (Intercept) 9.72013269 11.69401978 ## CO2 0.00732192 0.01280289

The 95% (default confidence level) confidence interval for the slope of CO2 is (0.0073, 0.0128). We are 95% confident that as CO2 concentration increases by 1 ppt the air temperature on average will increase by the amount between 0.0073 and 0.0128 degrees Celsius.

We now need to check if the fitted model meets the assumptions. From the scatterplot we can see the linearity assumption is satisfied. Because the data are a time series, we need to plot the residuals against the time (year) to check independence assumption. We also need to plot the residuals against fitted values ŷ (or x-values) to check constant variance assumptions. Let’s produce both plots: par(mfrow = c(1, 2)) plot(climate$year, imod$residuals, xlab = 'Time', ylab = 'Residual') abline(a = 0, b = 0)

2

plot(imod$fitted.values, imod$residuals, xlab = 'Fitted value', ylab = 'Residual') abline(a = 0, b = 0)

1980 1990 2000

− 0

.1 5

− 0

.0 5

0 .0

5 0

.1 5

Time

R e

si d

u a

l

14.1 14.3 14.5

− 0

.1 5

− 0

.0 5

0 .0

5 0

.1 5

Fitted value

R e

si d

u a

l

The two residual plots are almost identical. Both of them show no pattern and equal spread across the x-values, so the independence and constant variance assumptions are met. To check the Normal distribution assumption we may use histogram and Q-Q plot of residuals par(mfrow = c(1,2)) hist(imod$residuals, xlab = 'Residual', main = 'Histogram of Residuals') qqnorm(imod$residuals) qqline(imod$residuals)

Histogram of Residuals

Residual

F re

q u

e n

cy

−0.2 0.0 0.1 0.2

0 2

4 6

8

−2 −1 0 1 2

− 0

.1 5

− 0

.0 5

0 .0

5 0

.1 5

Normal Q−Q Plot

Theoretical Quantiles

S a

m p

le Q

u a

n til

e s

The histogram shows a potential right-skewness in the distribution. The Q-Q plot shows a slight departure from the straight line especially at two ends. The assumption is not seriously violated, but the regression

3

analysis should proceed with caution. Note that it is not easy to check the normality assumption with only 29 observations. More data should be collected for this purpose.

Now we want to predict the air temperature when the CO2 concentration level is at 355 ppt. First, we need to make sure that the 355 ppt is in the range of sampled x-values range(climate$CO2)

## [1] 338.67 384.84

Then, we create a new dataset for prediction new <- data.frame(CO2 = c(355))

Suppose we are interested in the mean temperature of all years with the CO2 concentration at 355 ppt. We therefore build a 95% confidence interval for the mean temperature: predict(imod, newdata = new, interval = 'confidence', level = 0.95)

## fit lwr upr ## 1 14.27923 14.2394 14.31906

The predicted mean temperature ŷν = 14.279 and its confidence interval is (14.239, 14.319). We are 95% confident that the mean air temperature of all years with CO2 concentration at 355 ppt is between 14.239 and 14.319 (degrees Celsius).

Suppose the CO2 concentration is estimated to be 355 ppt next year. What are the predicted temperature and the 95% prediction interval for that temperature for next year? predict(imod, newdata = new, interval = 'prediction', level = 0.95)

## fit lwr upr ## 1 14.27923 14.07329 14.48517

The predicted temperature for next year ŷν = 14.279, and its prediction interval is (14.073, 14.485). We are 95% confident that the temperature for next year is between 14.073 and 14.485 (degrees Celsius) if the CO2 concentration is at 355 ppt. Note that the prediction interval is wider than the confidence interval.

It is interesting to build a confidence band and a prediction band, and plot them to see how the predicted values and their uncertainties dynamically change with all possible x-values. To do so, we first need to create a new sequence of possible x-values within the range of sampled data: xx <- seq(min(climate$CO2), max(climate$CO2), length.out = 100)

and create a new dataset for prediction based on xx new.band <- data.frame(CO2 = xx)

Then, we build confidence interval and prediction interval for each value in new.band conf <- predict(imod, newdata = new.band, interval = 'confidence', level = .95) pred <- predict(imod, newdata = new.band, interval = 'prediction', level = .95)

To make the plot we first draw the data points and the regression line, and then add confidence and prediction bands to the plot plot(climate$CO2, climate$Temp, xlab = 'CO2 Concentration', ylab = 'Temperature') abline(imod) lines(xx, conf[, 'lwr'], lty = 2, col = 'red') lines(xx, conf[, 'upr'], lty = 2, col = 'red') lines(xx, pred[, 'lwr'], lty = 3, col = 'blue') lines(xx, pred[, 'upr'], lty = 3, col = 'blue')

4

340 350 360 370 380

1 4

.1 1

4 .3

1 4

.5

CO2 Concentration

Te m

p e

ra tu

re

As we can see, the prediction band is much wider than the confidence band, and it covers all the data points. Note that the prediction interval is NOT parallel, although it seems to be so. Both intervals actually get narrower as x-values approach their average, but wider as the x-values move away from it.

5