R Notebook

LAB REPORT 1

Name: GAN ZHI XUAN NO MATRIC: SD23053 SECTION: 02G

1. Extract data from the following excel file: gmp.txt.

file.choose()

[1] "C:\\Users\\User10\\Desktop\\Statistical Modelling And Simulation\\Lab Report 1\\gmp.txt"

gmp <- read.table("C:\\Users\\User10\\Desktop\\Statistical Modelling And Simulation\\Lab Report 1\\gmp.txt", header=TRUE)
gmp

2. Develop a multiple linear regression model by considering all variables.

model<-lm(y~x1+x2+x3+x4+x5+x6+x7+x8+x9+x10+x11, data=gmp)
model

Call:
lm(formula = y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + 
    x10 + x11, data = gmp)

Coefficients:
(Intercept)           x1           x2           x3           x4           x5           x6  
  17.339838    -0.075588    -0.069163     0.115117     1.494737     5.843495     0.317583  
         x7           x8           x9          x10          x11  
  -3.205390     0.180811    -0.397945    -0.005115     0.638483  

summary(model)

Call:
lm(formula = y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + 
    x10 + x11, data = gmp)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.3441 -1.6711 -0.4486  1.4906  5.2508 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept) 17.339838  30.355375   0.571   0.5749  
x1          -0.075588   0.056347  -1.341   0.1964  
x2          -0.069163   0.087791  -0.788   0.4411  
x3           0.115117   0.088113   1.306   0.2078  
x4           1.494737   3.101464   0.482   0.6357  
x5           5.843495   3.148438   1.856   0.0799 .
x6           0.317583   1.288967   0.246   0.8082  
x7          -3.205390   3.109185  -1.031   0.3162  
x8           0.180811   0.130301   1.388   0.1822  
x9          -0.397945   0.323456  -1.230   0.2344  
x10         -0.005115   0.005896  -0.868   0.3971  
x11          0.638483   3.021680   0.211   0.8350  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.227 on 18 degrees of freedom
  (因为不存在，2个观察量被删除了)
Multiple R-squared:  0.8355,    Adjusted R-squared:  0.7349 
F-statistic:  8.31 on 11 and 18 DF,  p-value: 5.231e-05

Part A

1. Construct a normal probability plot of the residuals. Does there seem to be any problems with the normality assumption?

model.stdres = rstandard(model)
qqnorm(model.stdres, ylab="Standardized Residuals",  xlab="Normal Scores",  main="Gallon") 
qqline(model.stdres)

Slight deviation appears at the tails, but most points lie close to the straight line, suggesting that residuals are approximately normal.

hist(residuals(model), col = "steelblue", main="Histogram of Residuals", xlab="Residuals")

The shape appears approximately bell-shaped. The residuals are approximately normally distributed.

2. Support your answer in (i) by using appropriate normality test.

shapiro.test(model.stdres)

    Shapiro-Wilk normality test

data:  model.stdres
W = 0.96009, p-value = 0.3113

W is a test statistics values.

\(H_{0}\): The data is normally distributed.

\(H_{1}\): The data is not normally distributed.

\(p-value=0.3113\)

Since \((p-value=0.3113)>(\alpha=0.05)\), do not reject \(H_{0}\)

At \(\alpha=0.05\), the data is normally distributed.

3. Construct and interpret a plot of residuals versus the predicted response.

yhat<- predict(model)   
epsilon <- residuals(model)
plot(cbind(yhat,epsilon), xlab="Predicted values", ylab="Residuals", main="Residuals vs Predicted")
abline(h = 0, lty = 2)

Interpretation:

-A horizontal band of residuals with no pattern. Therefore, linearity and homoscedasticity likely hold.

-The assumption of linearity appears satisfied. The residuals appear randomly distributed across the range of fitted values without any systematic curve or trend.

-The independence assumption appears reasonable. The residuals appear independent, showing no sequential structure or visible dependency.

-The scatter plot alone doesn’t clearly indicate normality.

-The equality of variances assumption appears satisfied. The spread of residuals appears quite uniform across the range of yhat and no funneling or widening pattern.

Part B

Then, detect any outliers occurs using Cook’s Distance method.

pairs(gmp)

round(cor(gmp),4)

          y      x1      x2 x3      x4      x5      x6      x7      x8      x9     x10     x11
y    1.0000 -0.8788 -0.8069 NA  0.3565  0.5948 -0.4870  0.7220 -0.7546 -0.7731 -0.8629 -0.7451
x1  -0.8788  1.0000  0.9452 NA -0.3302 -0.6316  0.6591 -0.7815  0.8552  0.8014  0.9457  0.8354
x2  -0.8069  0.9452  1.0000 NA -0.2921 -0.5170  0.7719 -0.6432  0.7974  0.7176  0.8834  0.7267
x3       NA      NA      NA  1      NA      NA      NA      NA      NA      NA      NA      NA
x4   0.3565 -0.3302 -0.2921 NA  1.0000  0.3737 -0.0493  0.4938 -0.2581 -0.3188 -0.2772 -0.3684
x5   0.5948 -0.6316 -0.5170 NA  0.3737  1.0000 -0.2054  0.8429 -0.5481 -0.4344 -0.5424 -0.7032
x6  -0.4870  0.6591  0.7719 NA -0.0493 -0.2054  1.0000 -0.3006  0.4252  0.3157  0.5206  0.4173
x7   0.7220 -0.7815 -0.6432 NA  0.4938  0.8429 -0.3006  1.0000 -0.6631 -0.6682 -0.7178 -0.8550
x8  -0.7546  0.8552  0.7974 NA -0.2581 -0.5481  0.4252 -0.6631  1.0000  0.8850  0.9476  0.6863
x9  -0.7731  0.8014  0.7176 NA -0.3188 -0.4344  0.3157 -0.6682  0.8850  1.0000  0.9015  0.6507
x10 -0.8629  0.9457  0.8834 NA -0.2772 -0.5424  0.5206 -0.7178  0.9476  0.9015  1.0000  0.7722
x11 -0.7451  0.8354  0.7267 NA -0.3684 -0.7032  0.4173 -0.8550  0.6863  0.6507  0.7722  1.0000

cooksd <- cooks.distance(model)

1. Plot the influential observation by Cook’s Distance and comments.

plot(cooksd, pch="*", cex=2, main="Influential Obs by Cooks distance") 
abline(h = 4*mean(cooksd, na.rm=T), col="red")  

text(x=1:length(cooksd)+1, y=cooksd, labels=ifelse(cooksd>4*mean(cooksd, na.rm=T),names(cooksd),""), col="red")

Interpretation:

There are two outliers which are 14 and 17. Observations 14 and 17 exceed the Cook’s distance cutoff, so they are influential. A point is considered influential if Cook’s D > (4×mean(D)).

2. Examine and list the point observation(s) which consider as outlier(s). Justify your answer.

influential <- as.numeric(names(cooksd)[(cooksd > 4*mean(cooksd, na.rm=T))])  
head(gmp[influential, ])

Interpretation:

The point observations are row 14 and 17.

Row 14 has very high y. Row 17 has very high x10. These points likely influence model estimates due to their extreme response or predictor values. Further examination may be needed to decide if they should be removed.

Part C

1. By considering only x1, x2, x3, x8, x9 and x10, construct the lack of fit test. Interpret and justify your answer.

sub_model<-lm(y~x1+x2+x3+x8+x9+x10, data=gmp)
summary(sub_model)

Call:
lm(formula = y ~ x1 + x2 + x3 + x8 + x9 + x10, data = gmp)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.7829 -1.6308 -0.2023  1.7894  6.2575 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) 31.891090  19.188065   1.662    0.110
x1          -0.051858   0.044919  -1.154    0.260
x2           0.001803   0.056564   0.032    0.975
x3           0.031761   0.070140   0.453    0.655
x8           0.129341   0.116709   1.108    0.279
x9          -0.206554   0.275463  -0.750    0.461
x10         -0.003947   0.004986  -0.792    0.437

Residual standard error: 3.206 on 23 degrees of freedom
  (因为不存在，2个观察量被删除了)
Multiple R-squared:  0.7924,    Adjusted R-squared:  0.7383 
F-statistic: 14.64 on 6 and 23 DF,  p-value: 7.75e-07

anova(model,sub_model)

Analysis of Variance Table

Model 1: y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10 + x11
Model 2: y ~ x1 + x2 + x3 + x8 + x9 + x10
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     18 187.40                           
2     23 236.43 -5   -49.024 0.9418  0.478

\(H_{0}\): Full model do not offers a statistically better fit than the reduced model.

\(H_{1}\): Full model offers a statistically better fit than the reduced model.

\(p-value=0.4780\)

Since \((p-value=0.4780)>(\alpha=0.05)\), do not reject \(H_{0}\)

At \(\alpha=0.05\), full model do not offers a statistically better fit than the reduced model.