Data

data <- as.data.frame(read_excel("D:/Ani/BSDS/Code/MDR/Practical Exam 1/data.xlsx"))
n <- nrow(data)
summary(data)
##        x1              x2               x3               x4        
##  Min.   : 44.7   Min.   : 2.700   Min.   : 4.600   Min.   : 1.450  
##  1st Qu.: 60.6   1st Qu.: 4.975   1st Qu.: 8.705   1st Qu.: 3.390  
##  Median : 73.2   Median : 7.080   Median :12.140   Median : 4.460  
##  Mean   : 76.3   Mean   : 7.323   Mean   :13.522   Mean   : 4.652  
##  3rd Qu.: 85.3   3rd Qu.: 9.115   3rd Qu.:16.920   3rd Qu.: 5.685  
##  Max.   :157.2   Max.   :18.080   Max.   :33.070   Max.   :11.960  
##        x5               x6               y         
##  Min.   : 4.350   Min.   : 4.050   Min.   : 41785  
##  1st Qu.: 7.875   1st Qu.: 7.975   1st Qu.: 59857  
##  Median :11.110   Median : 9.550   Median : 69177  
##  Mean   :12.002   Mean   :12.836   Mean   : 77756  
##  3rd Qu.:14.975   3rd Qu.:16.545   3rd Qu.: 92206  
##  Max.   :24.850   Max.   :43.370   Max.   :146345

43 observations, 6 predictors (x1-x6), response y.

Step 1: Correlations (check for collinearity)

round(cor(data), 3)
##       x1    x2    x3    x4    x5    x6     y
## x1 1.000 0.816 0.101 0.900 0.106 0.093 0.324
## x2 0.816 1.000 0.108 0.828 0.154 0.122 0.291
## x3 0.101 0.108 1.000 0.030 0.919 0.943 0.644
## x4 0.900 0.828 0.030 1.000 0.106 0.040 0.287
## x5 0.106 0.154 0.919 0.106 1.000 0.865 0.666
## x6 0.093 0.122 0.943 0.040 0.865 1.000 0.511
## y  0.324 0.291 0.644 0.287 0.666 0.511 1.000

x1, x2, x4 are pairwise correlated (r = 0.82-0.90). x3, x5, x6 are pairwise correlated (r = 0.86-0.94). These two blocks carry overlapping information, so using all six predictors together is expected to cause collinearity.

Step 2: Full model

full_model <- lm(y ~ x1 + x2 + x3 + x4 + x5 + x6, data = data)
summary(full_model)
## 
## Call:
## lm(formula = y ~ x1 + x2 + x3 + x4 + x5 + x6, data = data)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -32061  -9323   -922   8234  67823 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  12166.0    12526.9   0.971  0.33793   
## x1             183.7      310.1   0.593  0.55720   
## x2             111.5     1638.8   0.068  0.94612   
## x3            4353.5     1740.1   2.502  0.01704 * 
## x4            1116.8     3520.6   0.317  0.75290   
## x5            2016.6     1482.7   1.360  0.18226   
## x6           -2922.4     1067.3  -2.738  0.00955 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 17470 on 36 degrees of freedom
## Multiple R-squared:  0.5983, Adjusted R-squared:  0.5313 
## F-statistic: 8.936 on 6 and 36 DF,  p-value: 5.414e-06

The overall F-test is significant (p < 0.001), so the six predictors are jointly useful. But only x3 and x6 are individually significant. A significant F with mostly non-significant t-tests, on predictors already known to be correlated, is a sign that collinearity is inflating the standard errors of the redundant predictors.

Step 3: Variable selection

best_sub <- ols_step_best_subset(full_model)
best_sub$metrics[, c("n", "predictors", "adjr", "cp", "aic", "sbc")]
##   n        predictors      adjr        cp      aic      sbc
## 1 1                x5 0.4294405 10.914826 974.5088 979.7924
## 2 2             x1 x5 0.4832641  7.103534 971.1864 978.2312
## 3 3          x3 x4 x6 0.5445121  2.904106 966.6727 975.4787
## 4 4       x1 x3 x5 x6 0.5543582  3.133859 966.6161 977.1833
## 5 5    x1 x3 x4 x5 x6 0.5439506  5.004631 968.4620 980.7904
## 6 6 x1 x2 x3 x4 x5 x6 0.5313429  7.000000 970.4565 984.5461

Judging only by Adjusted R2, Mallow’s Cp, AIC and SBC: the model with x3, x4, x6 has the lowest AIC and SBC, a Cp value close to its own number of parameters (4), and an Adjusted R2 competitive with any larger model - so it is the preferred subset.

reduced_model <- lm(y ~ x3 + x4 + x6, data = data)
summary(reduced_model)
## 
## Call:
## lm(formula = y ~ x3 + x4 + x6, data = data)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -33249 -10851   -753   8634  67793 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    19644       9002   2.182  0.03518 *  
## x3              5968       1255   4.754 2.71e-05 ***
## x4              3463       1296   2.672  0.01095 *  
## x6             -3015       1042  -2.893  0.00622 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 17220 on 39 degrees of freedom
## Multiple R-squared:  0.577,  Adjusted R-squared:  0.5445 
## F-statistic: 17.74 on 3 and 39 DF,  p-value: 2.026e-07

Step 4: Non-linearity

plot(data$x3, resid(reduced_model), xlab = "x3", ylab = "Residuals")
abline(h = 0, lty = 2)

plot(data$x4, resid(reduced_model), xlab = "x4", ylab = "Residuals")
abline(h = 0, lty = 2)

plot(data$x6, resid(reduced_model), xlab = "x6", ylab = "Residuals")
abline(h = 0, lty = 2)

No curved pattern against any predictor, so a linear form is adequate.

Step 5: Outliers

r <- resid(reduced_model)
s <- summary(reduced_model)$sigma
plot(r, xlab = "Index", ylab = "Residuals")
abline(h = c(-2 * s, 2 * s), lty = 2)

r[abs(r) > 2 * s]
##       13 
## 67792.69
boxplot(r, main = "Residuals (raw y)")

Observation 13 has a residual well beyond 2 residual standard deviations, so its y-value is poorly predicted by its x’s. It does not look like a data error, so it is kept rather than removed.

Step 6: Normality of residuals

hist(resid(reduced_model), main = "Residuals (raw y)", xlab = "Residuals")

qqnorm(resid(reduced_model)); qqline(resid(reduced_model))

The histogram is right-skewed and the QQ-plot departs from the line in the tails, matching the skew already visible in y (median 69,177 vs mean 77,756, max 146,345).

Step 7: Heteroscedasticity, and the transformed final model

plot(fitted(reduced_model), resid(reduced_model), xlab = "Fitted", ylab = "Residuals")
abline(h = 0, lty = 2)

No strong funnel shape here, but the skewness and non-normal residuals above are reason enough to try a log transform of y:

final_model <- lm(log(y) ~ x3 + x4 + x6, data = data)
summary(final_model)
## 
## Call:
## lm(formula = log(y) ~ x3 + x4 + x6, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.43246 -0.15439  0.00942  0.10427  0.68346 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 10.50931    0.11018  95.379  < 2e-16 ***
## x3           0.07124    0.01537   4.636 3.92e-05 ***
## x4           0.03823    0.01587   2.409   0.0208 *  
## x6          -0.03414    0.01276  -2.676   0.0108 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2108 on 39 degrees of freedom
## Multiple R-squared:  0.5804, Adjusted R-squared:  0.5481 
## F-statistic: 17.98 on 3 and 39 DF,  p-value: 1.741e-07
confint(final_model)
##                    2.5 %      97.5 %
## (Intercept) 10.286437521 10.73217583
## x3           0.040156062  0.10231592
## x4           0.006136352  0.07032492
## x6          -0.059950615 -0.00833391
hist(resid(final_model), main = "Residuals (log y)", xlab = "Residuals")

qqnorm(resid(final_model)); qqline(resid(final_model))

plot(fitted(final_model), resid(final_model), xlab = "Fitted (log scale)", ylab = "Residuals")
abline(h = 0, lty = 2)

On the log scale the QQ-plot hugs the reference line, the histogram is roughly symmetric, and the residual-vs-fitted plot is an even band with no funnel. Both normality and constant variance look reasonable here.

Final model and interpretation

anova(final_model)
## Analysis of Variance Table
## 
## Response: log(y)
##           Df  Sum Sq Mean Sq F value   Pr(>F)    
## x3         1 1.84031 1.84031 41.4103 1.29e-07 ***
## x4         1 0.23862 0.23862  5.3694  0.02583 *  
## x6         1 0.31820 0.31820  7.1602  0.01084 *  
## Residuals 39 1.73319 0.04444                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model: log(y) = 10.509 + 0.0712 x3 + 0.0382 x4 - 0.0341 x6

Term Estimate Std. Error t p-value
Intercept 10.509 0.110 95.38 < 0.001
x3 0.0712 0.0154 4.64 < 0.001
x4 0.0382 0.0159 2.41 0.021
x6 -0.0341 0.0128 -2.68 0.011

Residual SE = 0.211 on 39 df. R-squared = 0.580. Adjusted R-squared = 0.548. F(3, 39) = 17.98, p < 0.001.

A caveat on x6’s sign. x6 is positively correlated with y on its own:

cor(data$x6, data$y)
## [1] 0.510605

but its partial coefficient in the final model is negative. This is not a contradiction, but it is a flag. x3 and x6 come from the same correlated block identified in Step 1 (r = 0.943 between them), and both remain in the final model together. When two correlated predictors are both in a model, it is common for one of them to pick up a partial coefficient with the opposite sign of its own simple correlation with y - the coefficient is measuring the part of x6 left over after netting out x3, not x6’s standalone relationship with y. So the “3.4% decrease” above should be read as the effect of x6 net of x3, not as evidence that raising x6 by itself would lower y. Because x3 and x6 still share a lot of variance in this model, that partial effect - and its sign - should be treated as less stable than the other estimates here, and not over-interpreted causally.

Conclusion

Starting from six predictors, two blocks of collinear variables were reduced to x3, x4 and x6 using best-subset selection judged on Adjusted R-squared, Cp, AIC and SBC. The relationship is linear in these predictors, one outlier (obs. 13) was identified and retained, and a log transform of y resolved the non-normality seen on the raw scale while keeping the residuals free of any strong non-constant-variance pattern. The final model explains about 58% of the variation in log(y) and all three predictors are statistically significant.