Statistical Modeling I
Practical in R – Output
Practical in R – Output
In this practical, we will work with the Stackloss dataset ([Link]). We will look at two
different models and their analysis.
The data are obtained in a production process of oxidizing ammonia. The variables of interest
are:
• Y : the stack loss, which is the percentage of the ingoing ammonia that escapes unab-
sorbed;
• X1 : the airflow
• X2 : the cooling water inlet temperature in degrees C;
• X3 : the acid concentration in percent.
i.i.d.
1. Fit Model 1: Yi = β0 + β1 x1,i + βx2,i + β3 x3,i + εi with εi ∼ N (0, σ 2 )
First of all, we load the data by using the usual command:
> data <- [Link](’[Link]’)
> attach(data)
Then we fit the model 1 by using the following commands:
> [Link] <- lm(y ~ x1 + x2 + x3)
> summary([Link])
Call:
lm(formula = y ~ x1 + x2 + x3)
Residuals:
Min 1Q Median 3Q Max
-0.72377 -0.17117 -0.04551 0.23614 0.56978
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.61414 8.90214 0.406 0.68982
x1 0.07156 0.01349 5.307 5.8e-05 ***
x2 0.12953 0.03680 3.520 0.00263 **
1
x3 -0.15212 0.15629 -0.973 0.34405
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.3243 on 17 degrees of freedom
Multiple R-squared: 0.9136,Adjusted R-squared: 0.8983
F-statistic: 59.9 on 3 and 17 DF, p-value: 3.016e-09
2. Check if there are any apparent problems with the residuals;
In order to check if there are problems with the residuals, we need to create them and
then checking the normality and linearity conditions:
> stdres1 <-rstandard([Link])
> fits1<-fitted([Link])
> plot(fits1,stdres1, main="Std res vs fits, stackloss")
> qqnorm(stdres1, main="Q-Q Plot, stackloss")
> qqline(stdres1)
In Figure 1.1, we show the standardized residuals versus the fitted values (left) and the
QQ plot (right). For the normality probability plot, we do not contradict the normal-
ity assumption, which is confirmed by the Shapiro-Wilk test (p-value equal to 0.6451).
Regarding the residuals versus the fitted values, the plot casts some doubt on the as-
sumption of a constant variance but not the linearity of the model
Std res vs fits, stackloss Q−Q Plot, stackloss
2
2
1
1
Sample Quantiles
stdres1
0
−1
−1
−2
−2
0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 −2 −1 0 1 2
fits1 Theoretical Quantiles
Figure 1.1: Plot of standardized residuals versus fitted values (left) and of the QQ plot for the
model with three explanatory variables (right).
2
3. Test the hypothesis regarding the overall regression by using the F test
Moving to the F test, we look at the last line of the summary in question 1), which is
F-statistic: 59.9 on 3 and 17 DF, p-value: 3.016e-09
Thus there is strong evidence against the null hypothesis (F = 59.9):
H0 : β1 = β2 = β3 = 0 versus H1 : one of the βi is different from zero
4. Test the hypothesis regarding the parameters βj for j = 0, 1, 2, 3 by using the t tests
Moving to the parameters, from the summary of the fitted models, we can see the t-
tests for each coefficient. In particular, β1 and β2 , the parameters related to x1 and
x2 are statistically significant, thus we reject the null hypothesis βj = 0 against the
alternative H1 : βj 6= 0. This conclusion is not confirmed for the intercept and the
parameter related to x3 , thus β0 and β3 are not statistically significant and we cannot
reject the null hypothesis. Obviously all these considerations are done when all the
other parameters are int the model.
i.i.d.
5. Fit Model 2: Yi = β0 + β1 x1,i + βx2,i + εi with εi ∼ N (0, σ 2 )
Moving to the second model, we use the usual command to fit the linear regression
model with two explanatory variables:
> [Link] <- lm(y ~ x1 + x2)
> summary([Link])
Call:
lm(formula = y ~ x1 + x2)
Residuals:
Min 1Q Median 3Q Max
-0.75290 -0.17505 0.01894 0.21156 0.56588
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -5.03588 0.51383 -9.801 1.22e-08 ***
x1 0.06712 0.01267 5.298 4.90e-05 ***
x2 0.12954 0.03675 3.525 0.00242 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.3239 on 18 degrees of freedom
Multiple R-squared: 0.9088,Adjusted R-squared: 0.8986
F-statistic: 89.64 on 2 and 18 DF, p-value: 4.382e-10
3
6. Check if there are any apparent problems with the residuals;
As done above, we need to define the standardized residuals in order to check the nor-
mality and linearity assumptions:
> stdres2 <-rstandard([Link])
> fits2<-fitted([Link])
> plot(fits2,stdres2, main="Std res vs fits, stackloss2")
> qqnorm(stdres2, main="Q-Q Plot, stackloss2")
> qqline(stdres2)
Figure 1.2 shows in the right panel that there is no problem with the normality assump-
tion, which is also confirmed from the Shapiro-Wilk test (p-value of 0.7321). Regarding
the standardized residuals versus the fitted values, the left panel shows some doubts on
the assumption of a constant variance but not regarding the linearity assumption.
Std res vs fits, stackloss2 Q−Q Plot, stackloss2
2
2
1
1
Sample Quantiles
0
0
stdres2
−1
−1
−2
−2
1.0 1.5 2.0 2.5 3.0 3.5 −2 −1 0 1 2
fits2 Theoretical Quantiles
Figure 1.2: Plot of standardized residuals versus fitted values (left) and of the QQ plot for the
model with two explanatory variables (right).
7. Test the hypothesis regarding the overall regression by using the F test
To test the overall regression, we use the F test from the summary of the linear regres-
sion and in this case the values of the F test is 89.64, thus the overall regression is highly
significant.
8. Test the hypothesis regarding the parameters βj for j = 0, 1, 2 by using the t tests
Moving to the parameters, from the summary of the linear regression, we can see that
all the three parameters are statistically significant and thus we can reject the null hy-
pothesis that the coefficients are equal to zero.
4
9. Which is the best model between Model 1 and Model 2 and why?
Looking at the adjusted R2 we can see that the Model 2 is better than Model 1 since
the adj(R2 ) = 89.86% that is high and comparable with the adjusted coefficient of
determination for the full model, which is adj(R2 ) = 89.83%