MATH1312 Regression Analysis Semester 1, 2021

Assignment3

1. (Data: “asphalt”) The data (n = 31) deals with pavement durability which contains measurements on the following variables: y= change in rut depth

x1= viscosity of rut depth

x2= % of asphalt in the surface course

x3= % of asphalt in the base course

x4= % of fines in the surface course

x5= % of voids in the surface course

x6= run indicator

x6 is a run indicator which separates the data into two different experimental runs:

1. Fit the full regression model with six predictors to the data set and use the ANOVA table to assess its overall fit.

2. Exhibit the fitted equation for y when the run indicator is 1 and when it is -1.

3. Use All Possible Subsets regression to select the best model based on their R2 values. Perform all diagnostic tests and check the adequacy of this model.

library(readr)
library(dplyr)
library(tidyverse)
library(car)

We have imported the necessary packages into R which will come handy in our further assignment.

#import the data
asphalt <- read.csv("asphalt.csv")
asphalt

Here we have imported the data into R and have a glimpse at our asphalt dataset.

1.a) Regression Model

asphaltmodel <- lm(y~visc+surf+base+fines+voids+run, asphalt)
summary(asphaltmodel)

Call:
lm(formula = y ~ visc + surf + base + fines + voids + run, data = asphalt)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.6781 -1.8309  0.1751  1.4858 11.1262 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -62.970450  36.118989  -1.743   0.0941 .  
visc          0.003071   0.008161   0.376   0.7100    
surf          7.498028   3.967155   1.890   0.0709 .  
base          6.225817   4.812723   1.294   0.2081    
fines         0.522211   1.174673   0.445   0.6606    
voids        -0.241275   1.684963  -0.143   0.8873    
run          -5.386297   0.985384  -5.466 1.28e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.94 on 24 degrees of freedom
Multiple R-squared:  0.7274,    Adjusted R-squared:  0.6592 
F-statistic: 10.67 on 6 and 24 DF,  p-value: 8.588e-06

Above we have fitted the regression model with all the six predictors available in the dataset by using Lm() Function and named it as asphaltmodel. Here, we can see that the R2 value is 0.7274, Adjusted R-square is 0.6592. To check the overall fit, we will use Anova function

anova(asphaltmodel)

Analysis of Variance Table

Response: y
          Df Sum Sq Mean Sq F value    Pr(>F)    
visc       1 424.93  424.93 27.3767 2.310e-05 ***
surf       1  16.13   16.13  1.0393    0.3182    
base       1  28.40   28.40  1.8295    0.1888    
fines      1  19.73   19.73  1.2711    0.2707    
voids      1  41.05   41.05  2.6450    0.1169    
run        1 463.77  463.77 29.8792 1.283e-05 ***
Residuals 24 372.51   15.52                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

After performing the Anova test, we found out that P value for visc and run variables are 0.001 which is less than 0.005. Values for all the other variables are greater then 0.005.

1.b) Fitting Equation.

prun <- filter(asphalt, run == "1")
prunmodel <- lm(y~visc+surf+base+fines+voids, prun)
summary(prunmodel)

Call:
lm(formula = y ~ visc + surf + base + fines + voids, data = prun)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.61544 -0.11773 -0.05099  0.13625  0.46611 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)  
(Intercept) -4.0419695  4.7947737  -0.843    0.421  
visc        -0.0024706  0.0008347  -2.960    0.016 *
surf         0.1755011  0.5703476   0.308    0.765  
base         0.3505746  0.8900850   0.394    0.703  
fines        0.1902141  0.2179108   0.873    0.405  
voids        0.2962432  0.2299413   1.288    0.230  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3703 on 9 degrees of freedom
Multiple R-squared:  0.6382,    Adjusted R-squared:  0.4372 
F-statistic: 3.175 on 5 and 9 DF,  p-value: 0.06313

nrun <- filter(asphalt, run == "-1")
nrunmodel <- lm(y~visc+surf+base+fines+voids, nrun)
summary(nrunmodel)

Call:
lm(formula = y ~ visc + surf + base + fines + voids, data = nrun)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.5155 -2.1402 -0.3626  1.8615  5.5877 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept) -147.0491    69.9536  -2.102   0.0619 .
visc          -1.7729     0.6298  -2.815   0.0183 *
surf          19.5080     7.5535   2.583   0.0273 *
base          10.4327     5.9136   1.764   0.1082  
fines          0.1704     1.4195   0.120   0.9068  
voids          3.6021     2.8706   1.255   0.2381  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.504 on 10 degrees of freedom
Multiple R-squared:  0.7356,    Adjusted R-squared:  0.6035 
F-statistic: 5.565 on 5 and 10 DF,  p-value: 0.01044

In the above step, we have subset our data into 2 category based on run variables, where values are 1 and -1. We need to find the fitted equation for Y variable with both the run values.

For that, we created prun data frame for run = 1 and nrun dataframe for run = -1. The fitted quation for prun is Y^ = -4.04 -0.02X1+ 0.17X2+ 0.35X3+ 0.19X4+ 0.29X5.
We can also see that the P-value is smaller then 0.05, therefore we can say that the model is significant and we can reject H0. Multiple R square is 0.638 and Adjusted R- square is observed as 0.437.

For *run = -1, we created a model based on dataset nrun. The fitted equation for nrun is Y^ = -147.04 -1.77X1 +19.50X2 +10.43X3 +0.17X4 +3.60X5**.

1.c) Use All Possible Subsets regression to select the best model based on their R2 values.

library(leaps)

package 㤼㸱leaps㤼㸲 was built under R version 4.0.4

r<-leaps::regsubsets(y~visc+surf+base+fines+voids+run, data=asphalt)
summary(r)

Subset selection object
Call: regsubsets.formula(y ~ visc + surf + base + fines + voids + run, 
    data = asphalt)
6 Variables  (and intercept)
      Forced in Forced out
visc      FALSE      FALSE
surf      FALSE      FALSE
base      FALSE      FALSE
fines     FALSE      FALSE
voids     FALSE      FALSE
run       FALSE      FALSE
1 subsets of each size up to 6
Selection Algorithm: exhaustive
         visc surf base fines voids run
1  ( 1 ) " "  " "  " "  " "   " "   "*"
2  ( 1 ) " "  "*"  " "  " "   " "   "*"
3  ( 1 ) " "  "*"  "*"  " "   " "   "*"
4  ( 1 ) " "  "*"  "*"  "*"   " "   "*"
5  ( 1 ) "*"  "*"  "*"  "*"   " "   "*"
6  ( 1 ) "*"  "*"  "*"  "*"   "*"   "*"

par(mfrow=c(1,1))
plot(r, scale="adjr2") 

Here, to check select the best possible regression model bases on R2, we have used leaps functions. Through the above graph, we can see that the best model based on R2 are surf, base and run.

1.c) All the diagnostic test to check the adequacy values.

To check the model adequacy, we are evaluating the assumptions of the model.

par(mfrow = c(2,2))
plot(asphaltmodel)

In the above plot we can see 4 Diagnostic plots known as Residuals plot, QQ plots, Scale-location plot and leverage plots.

In the 1st graph we can see that on X axis has predicted value known as Y^(y hat) and on y axis there are residuals and errors. Here we can the line is not flat and the points are looking as a bunch of clouds which means that linearity assumption is violated here.

In the next Graph (QQ plot), here y axis is ordered, observed and Standardized residuals and on X axis it has ordered theoretical residuals. In the graph we can see that residuals are truly normally distributed as in the beginning some were falling out of line but after that all the points are falling on a roughly straight line.

The next two shows that the regression is non-linear, non-constant variance.

To test further, if any assumptions are violated or not. We will perform some tests.

durbinWatsonTest(asphaltmodel)

 lag Autocorrelation D-W Statistic p-value
   1     -0.08199506       2.09818   0.968
 Alternative hypothesis: rho != 0

H0: Errors are uncorrelated H1: Errors are correlated

From the Durbin Waston test, we found out that p value is 0.984. This indictes negative/ No autocorrelation in the model. P value is grater then 0.05.

shapiro.test(asphaltmodel$residuals)

    Shapiro-Wilk normality test

data:  asphaltmodel$residuals
W = 0.93011, p-value = 0.04416

H0: Errors are normally distributed H1: Errors are not normally distributed

From the Shapiro-wilk normality test, after observing the P value of 0.015. We can say that the normality is violated here as p vale is greater then 0.05.

ncvTest(asphaltmodel)

Non-constant Variance Score Test 
Variance formula: ~ fitted.values 
Chisquare = 12.3959, Df = 1, p = 0.00043028

H0: Errors have a constant variance H1: Errors have a non-constant variance

Here, we can see that p value is less then 0.05. therefore, we can say that the constant variance assumptions is violated here.

2. (Data: “byssinosis”)

The data were collected from a group of workers in the cotton industry to assess the prevalence of the lung disease byssinosis among these workers. This disease is caused by long term exposure to particles of cotton, hemp, flax and jute working in this type of environment. It can result in asthma-like symptom which can lead to death among sufferers. The response variable y is binary and refer to number of workers suffering (response = yes) and not suffering (response = no) and the predictors are:

xl=dustiness of the workplace (1 = high, 2 = medium, 3 = low) x2=race ( 1 = European, 2 = other) x3=sex ( 1 = male, 2 = female) x4=smoking history (1 = smoker, 2 = nonsmoker) x5=length of employment in the cotton industry (1 = less than 10 years, 2 = between 10 and 20 years,3 = more than 20 years)

Notice that all five predictors are qualitative variables and the responses are entered in the event/trial format.

1. Fit a logistic regression model to the data set and discuss which of the predictors have a significant effect on the presence of byssinosis.
2. Discuss the adequacy of the logistic regression model.
3. (Needs to be done manually or by hand) From the final model you have selected in part (a), determine the probability that a person will suffer from byssinosis if given:

1. xl = 2, x2 = 2, x3 = 1, x4 = 2, x5 = 3,
2. xl = 1, x2 = 2, x3 = 2, x4 = 1, x5 = 3
3. xl = 2, x2 = 1, x3 = 1, x4 = 2, x5 = 2,
4. xl = 3, x2 = 1, x3 = 2, x4 = 2, x5 = 1.

#Q2
byssinosis <- read.csv("byssinosis.csv")
byssinosis

Here, we imported the byssinosis dataset that will be used in our Q2.

(a) Fit a logistic regression model to the data set and discuss which of the predictors have a significant effect on the presence of byssinosis.

#2.a)
glmmodel <- glm(cbind(BysYes,BysNo)~Dust+Race+Sex+Smoke+Employ,
                         family=binomial,byssinosis)
summary.glm(glmmodel)

Call:
glm(formula = cbind(BysYes, BysNo) ~ Dust + Race + Sex + Smoke + 
    Employ, family = binomial, data = byssinosis)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.4126  -0.7573  -0.2421   0.3688   1.9804  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -0.4852     0.6060  -0.801  0.42331    
Dust         -1.3751     0.1155 -11.901  < 2e-16 ***
Race          0.2463     0.2061   1.195  0.23203    
Sex          -0.2590     0.2116  -1.224  0.22095    
Smoke        -0.6292     0.1931  -3.259  0.00112 ** 
Employ        0.3856     0.1069   3.607  0.00031 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 322.527  on 64  degrees of freedom
Residual deviance:  69.509  on 59  degrees of freedom
AIC: 188.19

Number of Fisher Scoring iterations: 5

Here, we have used glm function to fit a logistic regression model and cbind function to bind the columns of Y variable, Family = binomial is used. In the above glm function, we can see that predictors Race, Smoke and Employ has significant impact on the presence of Byssinosis as their P-value is smaller then 0.05. Rest predictors has insignifacnt impat as they have higher P-value then 0.05

2.b) Discuss the adequacy of the logistic regression model.

pchisq(69.509, 59, lower.tail = FALSE)

[1] 0.1645647

To check the adequacy of the logistic regression model, we have used pchisq function by using the residual deviance and it’s **Degree of freedom*.

H0: the model fits the data well vs.  H1: the model does not fit the data well The chi-square test statistic of 69.509 with 59 degree of freedom gives a p-value of 0.1645, indicating that the null hypothesis is plausible, and we can conclude that logistic model is adequate.

2.c) determine the probability that a person will suffer from byssinosis if given:

i. xl = 2, x2 = 2, x3 = 1, x4 = 2, x5 = 3

1/(1+exp((-1.37*2)+(0.24*2)+(-0.25*1)+(-0.62*2)+(0.38*3)))

[1] 0.9315024

The probability that a person will suffer from byssinosis if xl = 2, x2 = 2, x3 = 1, x4 = 2, x5 = 3 is 0.9315

ii. xl = 1, x2 = 2, x3 = 2, x4 = 1, x5 = 3

1/(1+exp((-1.37*1)+(0.24*2)+(-0.25*2)+(-0.62*1)+(0.38*3)))

[1] 0.7047457

The probability that a person will suffer from byssinosis if xl = 1, x2 = 2, x3 = 2, x4 = 1, x5 = 3 is 0.7047

iii. xl = 2, x2 = 1, x3 = 1, x4 = 2, x5 = 2

1/(1+exp((-1.37*2)+(0.24*1)+(-0.25*1)+(-0.62*2)+(0.38*2)))

[1] 0.9619478

The probability that a person will suffer from byssinosis if xl = 2, x2 = 1, x3 = 1, x4 = 2, x5 = 2 is 0.9619

iv. xl = 3, x2 = 1, x3 = 2, x4 = 2, x5 = 1.

1/(1+exp((-1.37*3)+(0.24*1)+(-0.25*2)+(-0.62*2)+(0.38*1)))

[1] 0.994675

The probability that a person will suffer from byssinosis if xl = 3, x2 = 1, x3 = 2, x4 = 2, x5 = 1 is 0.9946