library(MASS)
#grab file and read in data
pressure <- read.csv(file.choose())
pressure
## x1 x2 x3 x4 y
## 1 2.14 10.00 0.34 1.000 28.9
## 2 4.14 10.00 0.34 1.000 31.0
## 3 8.15 10.00 0.34 1.000 26.4
## 4 2.14 10.00 0.34 0.246 27.2
## 5 4.14 10.00 0.34 0.379 26.1
## 6 8.15 10.00 0.34 0.474 23.2
## 7 2.14 10.00 0.34 0.141 19.7
## 8 4.14 10.00 0.34 0.234 22.1
## 9 8.15 10.00 0.34 0.311 22.8
## 10 2.14 10.00 0.34 0.076 29.2
## 11 4.14 10.00 0.34 0.132 23.6
## 12 8.15 10.00 0.34 0.184 23.6
## 13 2.14 2.63 0.34 0.679 24.2
## 14 4.14 2.63 0.34 0.804 22.1
## 15 8.15 2.63 0.34 0.890 20.9
## 16 2.14 2.63 0.34 0.514 17.6
## 17 4.14 2.63 0.34 0.672 15.7
## 18 8.15 2.63 0.34 0.801 15.8
## 19 2.14 2.63 0.34 0.346 14.0
## 20 4.14 2.63 0.34 0.506 17.1
## 21 8.15 2.63 0.34 0.669 18.3
## 22 2.14 2.63 0.34 1.000 33.8
## 23 4.14 2.63 0.34 1.000 31.7
## 24 8.15 2.63 0.34 1.000 28.1
## 25 5.60 1.25 0.34 0.848 18.1
## 26 5.60 1.25 0.34 0.737 16.5
## 27 5.60 1.25 0.34 0.651 15.4
## 28 5.60 1.25 0.34 0.554 15.0
## 29 4.30 2.63 0.34 0.748 19.1
## 30 4.30 2.63 0.34 0.682 16.2
## 31 4.30 2.63 0.34 0.524 16.3
## 32 4.30 2.63 0.34 0.472 15.8
## 33 4.30 2.63 0.34 0.398 15.4
## 34 5.60 10.10 0.25 0.789 19.2
## 35 5.60 10.10 0.25 0.677 8.4
## 36 5.60 10.10 0.25 0.590 15.0
## 37 5.60 10.10 0.25 0.523 12.0
## 38 5.60 10.10 0.34 0.789 21.9
## 39 5.60 10.10 0.34 0.677 21.3
## 40 5.60 10.10 0.34 0.590 21.6
## 41 5.60 10.10 0.34 0.523 19.8
## 42 4.30 10.10 0.34 0.741 21.6
## 43 4.30 10.10 0.34 0.617 17.3
## 44 4.30 10.10 0.34 0.524 20.0
## 45 4.30 10.10 0.34 0.457 18.6
## 46 2.40 10.10 0.34 0.615 22.1
## 47 2.40 10.10 0.34 0.473 14.7
## 48 2.40 10.10 0.34 0.381 15.8
## 49 2.40 10.10 0.34 0.320 13.2
## 50 5.60 10.10 0.55 0.789 30.8
## 51 5.60 10.10 0.55 0.677 27.5
## 52 5.60 10.10 0.55 0.590 25.2
## 53 5.60 10.10 0.55 0.523 22.8
## 54 2.14 112.00 0.34 0.680 41.7
## 55 4.14 112.00 0.34 0.803 33.7
## 56 8.15 112.00 0.34 0.889 29.7
## 57 2.14 112.00 0.34 0.514 41.8
## 58 4.14 112.00 0.34 0.672 37.1
## 59 8.15 112.00 0.34 0.801 40.1
## 60 2.14 112.00 0.34 0.306 42.7
## 61 4.14 112.00 0.34 0.506 48.6
## 62 8.15 112.00 0.34 0.668 42.4
#separate data into respective variables
x1 <- pressure$x1 #fluid velocity
x2 <- pressure$x2 #kinematic velocity
x3 <- pressure$x3 #mesh opening-cm
x4 <- pressure$x4 #gas velocity
y <- pressure$y #pressure drop
#see how each variable looks separately to see how they behave
plot(x1,y) #doesn't seem linear whatsoever
plot(x2,y) #could result in a linear regression line
plot(x3,y) #could result in a linear regression line
plot(x4,y) #seems like it might almost be a horizontal line
From just an initial glance at the provided data, x2 (kinematic velocity) and x3 (mesh opening-cm) appear to have what could be a linear correlation with the pressure drop.
#QUESTIONs 1,2,3,4
#trying our complete model with all interactions
model_f <- lm(y~(x1+x2+x3+x4)^2,data=pressure)
plot(model_f)
summary(model_f)
##
## Call:
## lm(formula = y ~ (x1 + x2 + x3 + x4)^2, data = pressure)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.4804 -3.0766 -0.6635 2.9625 12.2221
##
## Coefficients: (2 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 15.88376 23.17863 0.685 0.49616
## x1 0.18696 0.78447 0.238 0.81255
## x2 0.37921 0.06332 5.989 1.89e-07 ***
## x3 -11.99940 67.31148 -0.178 0.85919
## x4 -8.86442 35.62553 -0.249 0.80446
## x1:x2 0.01155 0.00869 1.329 0.18955
## x1:x3 NA NA NA NA
## x1:x4 -1.11525 1.14847 -0.971 0.33592
## x2:x3 NA NA NA NA
## x2:x4 -0.38547 0.11962 -3.222 0.00218 **
## x3:x4 72.85976 103.15353 0.706 0.48308
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.683 on 53 degrees of freedom
## Multiple R-squared: 0.7496, Adjusted R-squared: 0.7118
## F-statistic: 19.83 on 8 and 53 DF, p-value: 1.947e-13
There appear to be two interactions that are undefined (x1:x3 and x2:x3). The slope of the correlation between variables is not quite 0, but quite close to it.
#trying the model with no interactions
mod_one_fact <- lm(y~(x1+x2+x3+x4), data=pressure)
plot(mod_one_fact)
#QUESTIONs 5,7
summary(mod_one_fact)
##
## Call:
## lm(formula = y ~ (x1 + x2 + x3 + x4), data = pressure)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.9958 -3.3092 -0.2419 3.3924 10.5668
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.89453 4.32508 1.363 0.17828
## x1 -0.47790 0.34002 -1.406 0.16530
## x2 0.18271 0.01718 10.633 3.78e-15 ***
## x3 35.40284 11.09960 3.190 0.00232 **
## x4 5.84391 2.90978 2.008 0.04935 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.014 on 57 degrees of freedom
## Multiple R-squared: 0.6914, Adjusted R-squared: 0.6697
## F-statistic: 31.92 on 4 and 57 DF, p-value: 5.818e-14
#trying the model with only two variable interactions
#QUESTIONs 6,7
mod_two_fact <- lm(y~(x1:x2+x1:x3+x1:x4+x2:x3+x2:x4+x3:x4),data=pressure)
plot(mod_two_fact)
summary(mod_two_fact)
##
## Call:
## lm(formula = y ~ (x1:x2 + x1:x3 + x1:x4 + x2:x3 + x2:x4 + x3:x4),
## data = pressure)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.3385 -3.3185 -0.5334 3.1770 11.7028
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.987462 1.907001 6.810 7.71e-09 ***
## x1:x2 0.011500 0.008534 1.348 0.18333
## x1:x3 -0.194535 1.228698 -0.158 0.87478
## x1:x4 -0.755736 0.607862 -1.243 0.21904
## x2:x3 1.099938 0.178072 6.177 8.32e-08 ***
## x2:x4 -0.377924 0.115296 -3.278 0.00182 **
## x3:x4 41.624433 9.908353 4.201 9.82e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.607 on 55 degrees of freedom
## Multiple R-squared: 0.7485, Adjusted R-squared: 0.7211
## F-statistic: 27.29 on 6 and 55 DF, p-value: 7.867e-15
#using MASS library stepAIC function to get the ideal model
#QUESTION 8
#?stepAIC()
stepAIC(model_f,direction="backward") #function from MASS library
## Start: AIC=199.73
## y ~ (x1 + x2 + x3 + x4)^2
##
##
## Step: AIC=199.73
## y ~ x1 + x2 + x3 + x4 + x1:x2 + x1:x3 + x1:x4 + x2:x4 + x3:x4
##
##
## Step: AIC=199.73
## y ~ x1 + x2 + x3 + x4 + x1:x2 + x1:x4 + x2:x4 + x3:x4
##
## Df Sum of Sq RSS AIC
## - x3:x4 1 10.942 1173.4 198.31
## - x1:x4 1 20.682 1183.1 198.82
## <none> 1162.4 199.73
## - x1:x2 1 38.737 1201.2 199.76
## - x2:x4 1 227.751 1390.2 208.82
##
## Step: AIC=198.31
## y ~ x1 + x2 + x3 + x4 + x1:x2 + x1:x4 + x2:x4
##
## Df Sum of Sq RSS AIC
## - x1:x4 1 19.837 1193.2 197.35
## <none> 1173.4 198.31
## - x1:x2 1 38.709 1212.1 198.32
## - x2:x4 1 228.394 1401.8 207.34
## - x3 1 249.320 1422.7 208.26
##
## Step: AIC=197.35
## y ~ x1 + x2 + x3 + x4 + x1:x2 + x2:x4
##
## Df Sum of Sq RSS AIC
## - x1:x2 1 32.307 1225.5 197.01
## <none> 1193.2 197.35
## - x2:x4 1 220.026 1413.2 205.84
## - x3 1 252.209 1445.4 207.24
##
## Step: AIC=197.01
## y ~ x1 + x2 + x3 + x4 + x2:x4
##
## Df Sum of Sq RSS AIC
## - x1 1 11.262 1236.8 195.57
## <none> 1225.5 197.01
## - x2:x4 1 207.286 1432.8 204.69
## - x3 1 248.430 1473.9 206.45
##
## Step: AIC=195.57
## y ~ x2 + x3 + x4 + x2:x4
##
## Df Sum of Sq RSS AIC
## <none> 1236.8 195.57
## - x3 1 243.60 1480.4 204.72
## - x2:x4 1 245.68 1482.4 204.81
##
## Call:
## lm(formula = y ~ x2 + x3 + x4 + x2:x4, data = pressure)
##
## Coefficients:
## (Intercept) x2 x3 x4 x2:x4
## 1.5226 0.3806 34.5106 9.5247 -0.3047
#checking validity of ideal model
#QUESTIONS 9,10,11
best_mod <- lm(y~(x2+x3+x4+x2:x4), data=pressure)
plot(best_mod)
summary(best_mod)
##
## Call:
## lm(formula = y ~ (x2 + x3 + x4 + x2:x4), data = pressure)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.959 -3.358 -1.131 3.040 11.646
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.52261 4.03964 0.377 0.70763
## x2 0.38056 0.06084 6.255 5.47e-08 ***
## x3 34.51062 10.29961 3.351 0.00144 **
## x4 9.52471 2.96093 3.217 0.00214 **
## x2:x4 -0.30472 0.09056 -3.365 0.00137 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.658 on 57 degrees of freedom
## Multiple R-squared: 0.7336, Adjusted R-squared: 0.7149
## F-statistic: 39.24 on 4 and 57 DF, p-value: 9.297e-16
#getting new values for x2,x3,x4 for confidence and prediction intervals
newx2 <- seq(min(x2),max(x2),length.out=100)
newx3 <- seq(min(x3),max(x3),length.out=100)
newx4 <- seq(min(x4),max(x4),length.out=100)
#finding confidence and prediction intervals
#QUESTIONS 13,14,15,16
conf <- predict(best_mod, data.frame(x2=newx2,x3=newx3,x4=newx4), interval="confidence", level=.95)
pred <- predict(best_mod, data.frame(x2=newx2,x3=newx3,x4=newx4), interval="prediction", level=.95)
conf
## fit lwr upr
## 1 11.32090 7.422137 15.21966
## 2 11.90746 8.131027 15.68389
## 3 12.48765 8.829037 16.14627
## 4 13.06149 9.516152 16.60682
## 5 13.62896 10.192358 17.06556
## 6 14.19007 10.857646 17.52249
## 7 14.74481 11.512015 17.97761
## 8 15.29319 12.155464 18.43092
## 9 15.83521 12.788004 18.88242
## 10 16.37086 13.409650 19.33208
## 11 16.90016 14.020427 19.77989
## 12 17.42309 14.620369 20.22580
## 13 17.93965 15.209519 20.66978
## 14 18.44985 15.787933 21.11177
## 15 18.95369 16.355677 21.55171
## 16 19.45117 16.912829 21.98951
## 17 19.94228 17.459480 22.42508
## 18 20.42703 17.995732 22.85833
## 19 20.90542 18.521699 23.28913
## 20 21.37744 19.037506 23.71737
## 21 21.84310 19.543288 24.14291
## 22 22.30240 20.039189 24.56560
## 23 22.75533 20.525361 24.98530
## 24 23.20190 21.001958 25.40184
## 25 23.64211 21.469143 25.81507
## 26 24.07595 21.927076 26.22483
## 27 24.50343 22.375919 26.63095
## 28 24.92455 22.815829 27.03327
## 29 25.33930 23.246960 27.43165
## 30 25.74770 23.669458 27.82593
## 31 26.14972 24.083459 28.21599
## 32 26.54539 24.489091 28.60169
## 33 26.93469 24.886466 28.98292
## 34 27.31763 25.275684 29.35958
## 35 27.69421 25.656831 29.73158
## 36 28.06442 26.029977 30.09886
## 37 28.42827 26.395172 30.46136
## 38 28.78575 26.752453 30.81905
## 39 29.13688 27.101838 31.17191
## 40 29.48164 27.443325 31.51995
## 41 29.82003 27.776898 31.86317
## 42 30.15207 28.102522 32.20161
## 43 30.47774 28.420145 32.53533
## 44 30.79704 28.729698 32.86439
## 45 31.10999 29.031099 33.18887
## 46 31.41657 29.324250 33.50888
## 47 31.71678 29.609041 33.82453
## 48 32.01064 29.885349 34.13593
## 49 32.29813 30.153041 34.44322
## 50 32.57926 30.411979 34.74654
## 51 32.85402 30.662014 35.04603
## 52 33.12242 30.902995 35.34185
## 53 33.38446 31.134769 35.63416
## 54 33.64014 31.357181 35.92309
## 55 33.88945 31.570079 36.20882
## 56 34.13240 31.773313 36.49148
## 57 34.36898 31.966740 36.77123
## 58 34.59921 32.150222 37.04819
## 59 34.82307 32.323629 37.32250
## 60 35.04056 32.486843 37.59428
## 61 35.25170 32.639752 37.86364
## 62 35.45647 32.782257 38.13067
## 63 35.65487 32.914270 38.39547
## 64 35.84692 33.035713 38.65812
## 65 36.03260 33.146519 38.91867
## 66 36.21191 33.246633 39.17719
## 67 36.38487 33.336007 39.43373
## 68 36.55146 33.414606 39.68831
## 69 36.71169 33.482402 39.94097
## 70 36.86555 33.539373 40.19173
## 71 37.01305 33.585507 40.44059
## 72 37.15419 33.620798 40.68758
## 73 37.28896 33.645243 40.93268
## 74 37.41738 33.658846 41.17590
## 75 37.53942 33.661616 41.41723
## 76 37.65511 33.653562 41.65666
## 77 37.76443 33.634697 41.89417
## 78 37.86739 33.605039 42.12974
## 79 37.96399 33.564605 42.36337
## 80 38.05422 33.513412 42.59503
## 81 38.13809 33.451482 42.82470
## 82 38.21560 33.378835 43.05236
## 83 38.28674 33.295490 43.27799
## 84 38.35152 33.201470 43.50157
## 85 38.40994 33.096795 43.72308
## 86 38.46199 32.981486 43.94249
## 87 38.50768 32.855562 44.15980
## 88 38.54701 32.719045 44.37497
## 89 38.57997 32.571952 44.58799
## 90 38.60657 32.414304 44.79884
## 91 38.62681 32.246117 45.00750
## 92 38.64068 32.067410 45.21396
## 93 38.64820 31.878200 45.41819
## 94 38.64934 31.678502 45.62019
## 95 38.64413 31.468332 45.81993
## 96 38.63255 31.247705 46.01740
## 97 38.61461 31.016635 46.21259
## 98 38.59031 30.775135 46.40548
## 99 38.55964 30.523219 46.59606
## 100 38.52261 30.260898 46.78432
pred
## fit lwr upr
## 1 11.32090 1.211250 21.43054
## 2 11.90746 1.844354 21.97056
## 3 12.48765 2.468168 22.50714
## 4 13.06149 3.082809 23.04017
## 5 13.62896 3.688393 23.56952
## 6 14.19007 4.285035 24.09510
## 7 14.74481 4.872851 24.61677
## 8 15.29319 5.451955 25.13443
## 9 15.83521 6.022459 25.64796
## 10 16.37086 6.584476 26.15725
## 11 16.90016 7.138115 26.66220
## 12 17.42309 7.683484 27.16269
## 13 17.93965 8.220688 27.65861
## 14 18.44985 8.749830 28.14988
## 15 18.95369 9.271011 28.63637
## 16 19.45117 9.784328 29.11801
## 17 19.94228 10.289876 29.59468
## 18 20.42703 10.787744 30.06632
## 19 20.90542 11.278021 30.53281
## 20 21.37744 11.760792 30.99409
## 21 21.84310 12.236136 31.45006
## 22 22.30240 12.704129 31.90066
## 23 22.75533 13.164846 32.34581
## 24 23.20190 13.618354 32.78545
## 25 23.64211 14.064718 33.21950
## 26 24.07595 14.503998 33.64790
## 27 24.50343 14.936252 34.07061
## 28 24.92455 15.361532 34.48757
## 29 25.33930 15.779884 34.89873
## 30 25.74770 16.191353 35.30404
## 31 26.14972 16.595978 35.70347
## 32 26.54539 16.993793 36.09699
## 33 26.93469 17.384830 36.48455
## 34 27.31763 17.769114 36.86615
## 35 27.69421 18.146666 37.24174
## 36 28.06442 18.517504 37.61133
## 37 28.42827 18.881640 37.97489
## 38 28.78575 19.239082 38.33242
## 39 29.13688 19.589834 38.68392
## 40 29.48164 19.933896 39.02937
## 41 29.82003 20.271262 39.36880
## 42 30.15207 20.601921 39.70221
## 43 30.47774 20.925861 40.02961
## 44 30.79704 21.243063 40.35102
## 45 31.10999 21.553502 40.66647
## 46 31.41657 21.857152 40.97598
## 47 31.71678 22.153981 41.27959
## 48 32.01064 22.443953 41.57732
## 49 32.29813 22.727026 41.86923
## 50 32.57926 23.003156 42.15536
## 51 32.85402 23.272294 42.43575
## 52 33.12242 23.534385 42.71046
## 53 33.38446 23.789373 42.97955
## 54 33.64014 24.037195 43.24308
## 55 33.88945 24.277785 43.50111
## 56 34.13240 24.511073 43.75372
## 57 34.36898 24.736986 44.00098
## 58 34.59921 24.955445 44.24297
## 59 34.82307 25.166369 44.47976
## 60 35.04056 25.369673 44.71145
## 61 35.25170 25.565268 44.93812
## 62 35.45647 25.753063 45.15987
## 63 35.65487 25.932963 45.37678
## 64 35.84692 26.104868 45.58896
## 65 36.03260 26.268679 45.79651
## 66 36.21191 26.424293 45.99953
## 67 36.38487 26.571602 46.19813
## 68 36.55146 26.710499 46.39242
## 69 36.71169 26.840875 46.58250
## 70 36.86555 26.962617 46.76848
## 71 37.01305 27.075613 46.95049
## 72 37.15419 27.179748 47.12863
## 73 37.28896 27.274908 47.30302
## 74 37.41738 27.360977 47.47377
## 75 37.53942 27.437839 47.64101
## 76 37.65511 27.505380 47.80484
## 77 37.76443 27.563485 47.96538
## 78 37.86739 27.612038 48.12274
## 79 37.96399 27.650928 48.27705
## 80 38.05422 27.680042 48.42840
## 81 38.13809 27.699272 48.57691
## 82 38.21560 27.708509 48.72268
## 83 38.28674 27.707647 48.86583
## 84 38.35152 27.696585 49.00645
## 85 38.40994 27.675222 49.14465
## 86 38.46199 27.643462 49.28052
## 87 38.50768 27.601211 49.41415
## 88 38.54701 27.548381 49.54563
## 89 38.57997 27.484886 49.67506
## 90 38.60657 27.410644 49.80250
## 91 38.62681 27.325577 49.92804
## 92 38.64068 27.229613 50.05176
## 93 38.64820 27.122682 50.17371
## 94 38.64934 27.004720 50.29397
## 95 38.64413 26.875668 50.41259
## 96 38.63255 26.735469 50.52963
## 97 38.61461 26.584072 50.64515
## 98 38.59031 26.421431 50.75918
## 99 38.55964 26.247504 50.87177
## 100 38.52261 26.062252 50.98296