Statistic Basics and Linear Regression
This post explores some of the basic concepts of statistics. I mostly explore these concepts using linear regression. This is a reproducible example if you have R Studio just make sure you have installed the correct packages.
#http://r-statistics.co/Linear-Regression.html
#https://www.statmethods.net/stats/regression.html
#http://r-statistics.co/Statistical-Tests-in-R.html
#http://www.sthda.com/english/articles/40-regression-analysis/166-predict-in-r-model-predictions-and-confidence-intervals/
#Dr. Sager Utexas datasets
data <- read.table('AustinApartmentRents1.txt', header = TRUE)
summary(data)## Rent Area
## Min. : 399.0 Min. : 474
## 1st Qu.: 470.0 1st Qu.: 666
## Median : 535.0 Median : 755
## Mean : 572.3 Mean : 816
## 3rd Qu.: 638.8 3rd Qu.: 925
## Max. :1050.0 Max. :1864cor(data$Rent, data$Area)## [1] 0.8740597model <- lm(Rent ~ Area, data = data)
summary(model)##
## Call:
## lm(formula = Rent ~ Area, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -154.659 -50.882 8.189 54.874 148.207
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 160.18706 31.36081 5.108 3.8e-06 ***
## Area 0.50497 0.03685 13.702 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 68.86 on 58 degrees of freedom
## Multiple R-squared: 0.764, Adjusted R-squared: 0.7599
## F-statistic: 187.7 on 1 and 58 DF, p-value: < 2.2e-16test.Areas <- data.frame(Area = c (500,1000))
predict(model, newdata = test.Areas)## 1 2
## 412.6713 665.1556data1 <- read.table('AustinApartmentRents2.txt', header = TRUE)#A convient tool to see a lot of the initial data exploration
#https://towardsdatascience.com/simple-fast-exploratory-data-analysis-in-r-with-dataexplorer-package-e055348d9619
library(DataExplorer)## Warning: package 'DataExplorer' was built under R version 4.0.2plot_str(data1)
plot_missing(data1)
plot_histogram(data1)
plot_density(data1)
plot_correlation(data1)
plot_bar(data1)
#create_report(data1) #This creates an HTML report of all the above information and more#Confidence intervals around indivudual values
pred.int <- predict(model, interval = 'prediction')## Warning in predict.lm(model, interval = "prediction"): predictions on current data refer to _future_ responses#Confidence intervals around means
pred.conf <- predict(model, interval = 'confidence')
cbind(data,pred.int,pred.conf)## Rent Area fit lwr upr fit lwr upr
## 1 519 725 526.2893 387.1539 665.4247 526.2893 507.2700 545.3085
## 2 765 995 662.6308 523.0320 802.2296 662.6308 640.4747 684.7868
## 3 475 481 403.0769 261.9228 544.2310 403.0769 372.6213 433.5326
## 4 575 925 627.2830 488.0776 766.4884 627.2830 607.7583 646.8077
## 5 415 600 463.1682 323.2840 603.0524 463.1682 439.2801 487.0563
## 6 530 668 497.5061 358.1044 636.9078 497.5061 476.6278 518.3843
## 7 580 725 526.2893 387.1539 665.4247 526.2893 507.2700 545.3085
## 8 995 1421 877.7474 731.7844 1023.7104 877.7474 829.7031 925.7917
## 9 565 672 499.5259 360.1470 638.9048 499.5259 478.8005 520.2514
## 10 620 1025 677.7799 537.9544 817.6053 677.7799 654.2379 701.3218
## 11 450 781 554.5675 415.5703 693.5648 554.5675 536.5869 572.5481
## 12 520 800 564.1619 425.1837 703.1402 564.1619 546.3289 581.9950
## 13 495 870 599.5097 460.4795 738.5399 599.5097 581.2764 617.7431
## 14 420 700 513.6651 374.4284 652.9017 513.6651 493.9190 533.4112
## 15 575 800 564.1619 425.1837 703.1402 564.1619 546.3289 581.9950
## 16 425 620 473.2676 333.5438 612.9913 473.2676 450.3375 496.1977
## 17 770 1040 685.3544 545.4026 825.3061 685.3544 661.0735 709.6352
## 18 445 520 422.7707 282.0919 563.4495 422.7707 394.5998 450.9416
## 19 510 880 604.5594 465.5062 743.6127 604.5594 586.1509 622.9679
## 20 635 832 580.3209 441.3427 719.2991 580.3209 562.4884 598.1535
## 21 470 545 435.3949 294.9906 575.7993 435.3949 408.6285 462.1614
## 22 700 921 625.2631 486.0744 764.4518 625.2631 605.8580 644.6682
## 23 450 577 451.5539 311.4663 591.6416 451.5539 426.5018 476.6060
## 24 785 1080 705.5531 565.2224 845.8838 705.5531 679.1757 731.9306
## 25 485 710 518.7147 379.5215 657.9080 518.7147 499.2771 538.1524
## 26 415 605 465.6930 325.8504 605.5357 465.6930 442.0494 489.3367
## 27 399 680 503.5657 364.2305 642.9008 503.5657 483.1366 523.9948
## 28 585 730 528.8141 389.6960 667.9322 528.8141 509.9220 547.7063
## 29 525 687 507.1005 367.8016 646.3994 507.1005 486.9201 527.2809
## 30 495 703 515.1800 375.9568 654.4032 515.1800 495.5288 534.8311
## 31 505 672 499.5259 360.1470 638.9048 499.5259 478.8005 520.2514
## 32 445 660 493.4663 354.0171 632.9155 493.4663 472.2734 514.6593
## 33 565 755 541.4383 402.3922 680.4845 541.4383 523.0835 559.7931
## 34 650 810 569.2116 430.2377 708.1855 569.2116 551.4123 587.0109
## 35 515 611 468.7229 328.9288 608.5169 468.7229 445.3683 492.0774
## 36 470 705 516.1899 376.9755 655.4044 516.1899 496.6009 535.7789
## 37 470 564 444.9893 304.7778 585.2009 444.9893 419.2531 470.7255
## 38 700 1250 791.3978 648.7851 934.0105 791.3978 754.7720 828.0235
## 39 455 512 418.7310 277.9593 559.5026 418.7310 390.1001 447.3618
## 40 550 630 478.3173 338.6680 617.9666 478.3173 455.8452 500.7893
## 41 625 850 589.4103 450.4146 728.4061 589.4103 571.4413 607.3794
## 42 745 1156 743.9307 602.7129 885.1486 743.9307 713.1810 774.6805
## 43 540 932 630.8178 491.5816 770.0540 630.8178 611.0749 650.5607
## 44 650 755 541.4383 402.3922 680.4845 541.4383 523.0835 559.7931
## 45 595 1093 712.1177 571.6507 852.5847 712.1177 685.0246 739.2108
## 46 470 751 539.4185 400.3624 678.4745 539.4185 520.9890 557.8479
## 47 480 608 467.2080 327.3898 607.0261 467.2080 443.7095 490.7064
## 48 460 900 614.6588 475.5477 753.7699 614.6588 595.8181 633.4994
## 49 600 860 594.4600 455.4490 733.4711 594.4600 576.3734 612.5467
## 50 575 925 627.2830 488.0776 766.4884 627.2830 607.7583 646.8077
## 51 659 944 636.8774 497.5840 776.1708 636.8774 616.7351 657.0197
## 52 650 940 634.8575 495.5838 774.1312 634.8575 614.8518 654.8632
## 53 750 1048 689.3941 549.3715 829.4168 689.3941 664.7079 714.0803
## 54 455 474 399.5422 258.2967 540.7876 399.5422 368.6660 430.4184
## 55 430 700 513.6651 374.4284 652.9017 513.6651 493.9190 533.4112
## 56 605 921 625.2631 486.0744 764.4518 625.2631 605.8580 644.6682
## 57 929 1229 780.7934 638.5205 923.0664 780.7934 745.5138 816.0731
## 58 695 896 612.6389 473.5406 751.7372 612.6389 593.8932 631.3846
## 59 455 630 478.3173 338.6680 617.9666 478.3173 455.8452 500.7893
## 60 1050 1864 1101.4485 942.4198 1260.4772 1101.4485 1022.1188 1180.7782# 2. Regression line + confidence intervals
library("ggplot2")## Warning: package 'ggplot2' was built under R version 4.0.2mydata <- cbind(data, pred.int)
p <- ggplot(mydata, aes(Area, Rent)) +
geom_point() +
stat_smooth(method = lm)
# 3. Add prediction intervals
p + geom_line(aes(y = lwr), color = "red", linetype = "dashed")+
geom_line(aes(y = upr), color = "red", linetype = "dashed")
#T Test for samples
library(dplyr)
sample1 <- sample_n(data,40)
model1 <- lm(Rent ~ Area, data = data)
p1 <- predict(model1, interval = 'confidence', level = 0.95)
summary(p1)## fit lwr upr
## Min. : 399.5 Min. : 368.7 Min. : 430.4
## 1st Qu.: 496.5 1st Qu.: 475.5 1st Qu.: 517.5
## Median : 541.4 Median : 523.1 Median : 559.8
## Mean : 572.3 Mean : 548.8 Mean : 595.7
## 3rd Qu.: 627.3 3rd Qu.: 607.8 3rd Qu.: 646.8
## Max. :1101.4 Max. :1022.1 Max. :1180.8t.test(p1, mu = 550)##
## One Sample t-test
##
## data: p1
## t = 2.4118, df = 179, p-value = 0.01688
## alternative hypothesis: true mean is not equal to 550
## 95 percent confidence interval:
## 554.0485 590.4849
## sample estimates:
## mean of x
## 572.2667#MultiVariable Linear Regression
summary(data1)## Rent Area Bedrooms Bathrooms
## Min. : 399.0 Min. : 474 Min. :1.000 Min. :1.00
## 1st Qu.: 470.0 1st Qu.: 666 1st Qu.:1.000 1st Qu.:1.00
## Median : 535.0 Median : 755 Median :1.000 Median :1.00
## Mean : 572.3 Mean : 816 Mean :1.517 Mean :1.25
## 3rd Qu.: 638.8 3rd Qu.: 925 3rd Qu.:2.000 3rd Qu.:1.25
## Max. :1050.0 Max. :1864 Max. :5.000 Max. :2.00
## Security Parking Distance Shuttle
## Min. :0.0000 Min. :0.0000 Min. : 1.100 Min. :0.0000
## 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.: 5.000 1st Qu.:1.0000
## Median :0.0000 Median :0.0000 Median : 6.000 Median :1.0000
## Mean :0.1667 Mean :0.1333 Mean : 5.935 Mean :0.8667
## 3rd Qu.:0.0000 3rd Qu.:0.0000 3rd Qu.: 7.000 3rd Qu.:1.0000
## Max. :1.0000 Max. :1.0000 Max. :10.500 Max. :1.0000
## Age
## Min. : 1.00
## 1st Qu.:10.00
## Median :16.50
## Mean :16.33
## 3rd Qu.:22.25
## Max. :32.00model2 <- lm(Rent ~ Area + Bathrooms, data = data1)
summary(model2)##
## Call:
## lm(formula = Rent ~ Area + Bathrooms, data = data1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -152.02 -45.45 10.38 39.91 129.28
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 143.66927 29.51345 4.868 9.31e-06 ***
## Area 0.38746 0.04982 7.777 1.61e-10 ***
## Bathrooms 89.92902 27.75071 3.241 0.00199 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 63.83 on 57 degrees of freedom
## Multiple R-squared: 0.8007, Adjusted R-squared: 0.7937
## F-statistic: 114.5 on 2 and 57 DF, p-value: < 2.2e-16test.bathrooms <- data.frame(Area = c(500,1000), Bathrooms = c(1,2))
p2a <- predict(model2, newdata = test.bathrooms, interval = 'confidence')
p2 <- as.data.frame(predict(model2, interval = 'confidence', level = 0.95))
cbind(p2, data1)## fit lwr upr Rent Area Bedrooms Bathrooms Security Parking
## 1 514.5062 495.4259 533.5866 519 725 1 1 0 0
## 2 709.0493 673.7669 744.3316 765 995 2 2 0 0
## 3 419.9662 389.8581 450.0743 475 481 1 1 0 0
## 4 681.9271 643.6132 720.2411 575 925 2 2 0 1
## 5 466.0738 443.8498 488.2979 415 600 1 1 0 0
## 6 492.4211 472.8074 512.0348 530 668 1 1 0 0
## 7 514.5062 495.4259 533.5866 580 725 1 1 0 0
## 8 874.1069 829.4987 918.7151 995 1421 2 2 0 1
## 9 493.9709 474.4481 513.4937 565 672 1 1 0 0
## 10 630.7440 594.3939 667.0941 620 1025 2 1 1 0
## 11 536.2040 516.0356 556.3724 450 781 1 1 1 0
## 12 543.5657 522.6986 564.4328 520 800 2 1 0 0
## 13 570.6878 546.1304 595.2452 495 870 2 1 0 0
## 14 504.8198 485.7109 523.9287 420 700 1 1 0 0
## 15 543.5657 522.6986 564.4328 575 800 1 1 0 0
## 16 473.8230 452.5572 495.0888 425 620 1 1 0 0
## 17 726.4849 692.5303 760.4395 770 1040 2 2 0 1
## 18 435.0771 407.8700 462.2843 445 520 1 1 0 0
## 19 574.5624 549.3637 599.7611 510 880 2 1 0 1
## 20 555.9644 533.6045 578.3243 635 832 1 1 0 0
## 21 444.7636 419.2769 470.2502 470 545 1 1 0 0
## 22 680.3773 641.8590 718.8955 700 921 2 2 0 0
## 23 457.1623 433.6744 480.6501 450 577 1 1 0 0
## 24 741.9833 708.7412 775.2254 785 1080 2 2 0 0
## 25 508.6944 489.6360 527.7527 485 710 1 1 0 0
## 26 468.0111 446.0397 489.9826 415 605 1 1 0 0
## 27 497.0706 477.7062 516.4349 399 680 1 1 0 1
## 28 516.4435 497.3298 535.5573 585 730 2 1 0 0
## 29 499.7828 480.5310 519.0346 525 687 1 1 0 0
## 30 505.9821 486.8939 525.0704 495 703 1 1 0 0
## 31 493.9709 474.4481 513.4937 505 672 1 1 1 0
## 32 489.3214 469.5030 509.1398 445 660 1 1 0 0
## 33 526.1300 506.6576 545.6024 565 755 2 1 0 0
## 34 547.4403 526.1468 568.7337 650 810 2 1 0 0
## 35 470.3359 448.6563 492.0154 515 611 1 1 0 0
## 36 506.7571 487.6799 525.8342 470 705 1 1 0 0
## 37 452.1253 427.8562 476.3944 470 564 1 1 0 0
## 38 807.8514 772.3992 843.3035 700 1250 3 2 0 1
## 39 431.9774 404.1949 459.7599 455 512 1 1 1 0
## 40 477.6976 456.8558 498.5395 550 630 1 1 0 0
## 41 562.9387 539.5887 586.2886 625 850 2 1 1 0
## 42 771.4302 738.2367 804.6237 745 1156 3 2 0 0
## 43 594.7103 565.8488 623.5718 540 932 2 1 0 0
## 44 526.1300 506.6576 545.6024 650 755 1 1 1 1
## 45 747.0203 713.9093 780.1312 595 1093 2 2 1 0
## 46 524.5802 505.1862 543.9741 470 751 1 1 1 0
## 47 469.1735 447.3496 490.9974 480 608 1 1 0 0
## 48 582.3116 555.7642 608.8590 460 900 1 1 0 1
## 49 566.8132 542.8728 590.7537 600 860 2 1 0 0
## 50 591.9981 563.6574 620.3388 575 925 2 1 0 0
## 51 689.2888 651.9023 726.6754 659 944 2 2 0 0
## 52 687.7390 650.1632 725.3148 650 940 2 2 0 0
## 53 729.5846 695.8090 763.3602 750 1048 2 2 0 0
## 54 417.2540 386.6020 447.9060 455 474 1 1 0 0
## 55 504.8198 485.7109 523.9287 430 700 1 1 0 0
## 56 590.4483 562.4017 618.4949 605 921 1 1 0 0
## 57 799.7147 764.9734 834.4561 929 1229 2 2 1 0
## 58 670.6908 630.8290 710.5526 695 896 2 2 0 0
## 59 477.6976 456.8558 498.5395 455 630 1 1 1 0
## 60 1045.7513 964.5360 1126.9667 1050 1864 5 2 0 0
## Distance Shuttle Age
## 1 10.5 1 9
## 2 6.5 1 17
## 3 6.5 1 17
## 4 4.0 1 9
## 5 5.0 1 30
## 6 6.5 1 19
## 7 7.0 1 17
## 8 6.5 1 16
## 9 7.0 1 17
## 10 5.0 1 3
## 11 5.5 1 3
## 12 6.0 1 20
## 13 5.0 1 27
## 14 6.0 1 22
## 15 7.0 1 10
## 16 8.0 0 27
## 17 6.5 1 16
## 18 3.0 1 12
## 19 7.0 0 25
## 20 6.0 1 13
## 21 6.5 1 9
## 22 3.0 1 26
## 23 8.0 1 18
## 24 5.0 1 10
## 25 6.0 1 25
## 26 6.0 1 22
## 27 7.0 0 25
## 28 6.5 1 19
## 29 7.0 1 15
## 30 6.5 1 14
## 31 6.5 1 9
## 32 6.0 1 25
## 33 3.0 1 12
## 34 2.0 1 32
## 35 6.5 1 17
## 36 7.5 0 13
## 37 5.0 1 10
## 38 4.0 1 9
## 39 10.0 0 10
## 40 2.0 1 32
## 41 7.0 1 1
## 42 7.5 0 13
## 43 6.0 1 22
## 44 1.1 1 26
## 45 5.5 1 3
## 46 5.0 1 3
## 47 6.0 1 15
## 48 4.0 1 9
## 49 6.0 1 25
## 50 6.0 1 23
## 51 7.0 1 25
## 52 8.0 0 27
## 53 7.0 1 3
## 54 5.0 1 10
## 55 6.0 1 20
## 56 7.5 0 13
## 57 5.0 1 11
## 58 6.5 1 19
## 59 5.5 1 9
## 60 6.0 1 22p2a## fit lwr upr
## 1 427.3279 398.6624 455.9935
## 2 710.9866 675.8776 746.0955model3 <- lm(Rent ~ Area + Bathrooms + Security + Parking + Distance, data = data1)
summary(model3)##
## Call:
## lm(formula = Rent ~ Area + Bathrooms + Security + Parking + Distance,
## data = data1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -163.199 -37.278 4.548 38.345 149.276
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 167.45045 43.12346 3.883 0.000283 ***
## Area 0.39647 0.05084 7.798 2.09e-10 ***
## Bathrooms 92.60040 28.01804 3.305 0.001691 **
## Security -0.67875 22.42818 -0.030 0.975969
## Parking -38.26531 25.89638 -1.478 0.145316
## Distance -4.92937 5.05014 -0.976 0.333374
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 64.02 on 54 degrees of freedom
## Multiple R-squared: 0.81, Adjusted R-squared: 0.7925
## F-statistic: 46.06 on 5 and 54 DF, p-value: < 2.2e-16model4 <- lm(Rent ~ Distance + Parking + Security + Bathrooms + Area, data = data1)
summary(model4)##
## Call:
## lm(formula = Rent ~ Distance + Parking + Security + Bathrooms +
## Area, data = data1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -163.199 -37.278 4.548 38.345 149.276
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 167.45045 43.12346 3.883 0.000283 ***
## Distance -4.92937 5.05014 -0.976 0.333374
## Parking -38.26531 25.89638 -1.478 0.145316
## Security -0.67875 22.42818 -0.030 0.975969
## Bathrooms 92.60040 28.01804 3.305 0.001691 **
## Area 0.39647 0.05084 7.798 2.09e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 64.02 on 54 degrees of freedom
## Multiple R-squared: 0.81, Adjusted R-squared: 0.7925
## F-statistic: 46.06 on 5 and 54 DF, p-value: < 2.2e-16#Multicollinearity
#http://www.sthda.com/english/articles/39-regression-model-diagnostics/160-multicollinearity-essentials-and-vif-in-r/
library(caret)## Warning: package 'caret' was built under R version 4.0.2set.seed(1234)
y <- runif(50,min=0, max =100)
x1 <- runif(50, min = 0, max = 100)
x2 <- runif(50, min = 0, max = 100)
z1 <- x1+x2
z2 <- x1 + x2 + 0.005*runif(50,min=0, max =100)
list <- cbind(y,x1,x2,z1,z2)
list## y x1 x2 z1 z2
## [1,] 11.3703411 7.377988 3.545673 10.92366 11.14318
## [2,] 62.2299405 30.968660 56.507611 87.47627 87.59057
## [3,] 60.9274733 71.727174 28.025778 99.75295 99.79403
## [4,] 62.3379442 50.454591 20.419632 70.87422 71.29936
## [5,] 86.0915384 15.299896 13.373890 28.67379 28.79112
## [6,] 64.0310605 50.393349 32.568192 82.96154 83.45562
## [7,] 0.9495756 49.396092 15.506197 64.90229 65.20324
## [8,] 23.2550506 75.120020 12.996214 88.11623 88.61560
## [9,] 66.6083758 17.464982 43.553106 61.01809 61.20589
## [10,] 51.4251141 84.839241 3.864265 88.70351 88.98107
## [11,] 69.3591292 86.483383 71.330156 157.81354 158.02826
## [12,] 54.4974836 4.185728 10.076904 14.26263 14.55057
## [13,] 28.2733584 31.718216 95.030494 126.74871 126.96496
## [14,] 92.3433484 1.374994 12.181776 13.55677 13.66919
## [15,] 29.2315840 23.902573 21.965662 45.86823 45.91073
## [16,] 83.7295628 70.649462 91.308777 161.95824 162.27689
## [17,] 28.6223285 30.809476 94.585312 125.39479 125.61030
## [18,] 26.6820780 50.854757 27.915622 78.77038 78.80674
## [19,] 18.6722790 5.164662 12.347109 17.51177 17.91297
## [20,] 23.2225911 56.456984 79.716046 136.17303 136.33567
## [21,] 31.6612455 12.148019 74.427722 86.57574 86.95438
## [22,] 30.2693371 89.283638 91.597422 180.88106 181.17320
## [23,] 15.9046003 1.462726 99.459825 100.92255 101.27697
## [24,] 3.9995918 78.312110 94.236072 172.54818 172.76167
## [25,] 21.8799541 8.996133 48.613541 57.60967 57.78146
## [26,] 81.0598552 51.918998 28.345954 80.26495 80.64451
## [27,] 52.5697547 38.426669 25.154570 63.58124 63.79325
## [28,] 91.4658166 7.005250 50.325517 57.33077 57.61121
## [29,] 83.1345047 32.064442 49.696617 81.76106 81.81913
## [30,] 4.5770263 66.849540 31.844581 98.69412 98.84563
## [31,] 45.6091482 92.640048 96.222283 188.86233 189.10173
## [32,] 26.5186672 47.190972 63.409937 110.60091 110.77332
## [33,] 30.4672203 14.261534 12.743340 27.00487 27.30523
## [34,] 50.7306870 54.426976 42.304699 96.73167 96.76972
## [35,] 18.1096208 19.617465 91.431691 111.04916 111.52715
## [36,] 75.9670635 89.858049 46.779233 136.63728 136.64839
## [37,] 20.1248038 38.949978 90.816915 129.76689 130.18775
## [38,] 25.8809819 31.087078 59.774328 90.86141 91.17763
## [39,] 99.2150418 16.002866 63.174282 79.17715 79.33219
## [40,] 80.7352340 89.618585 86.915832 176.53442 176.90570
## [41,] 55.3333591 16.639378 50.274982 66.91436 67.23382
## [42,] 64.6406094 90.042460 98.363511 188.40597 188.90223
## [43,] 31.1824307 13.407820 32.438603 45.84642 45.91056
## [44,] 62.1819198 13.161413 48.137495 61.29891 61.74053
## [45,] 32.9770176 10.528750 35.698708 46.22746 46.63250
## [46,] 50.1997473 51.158358 62.747768 113.90613 114.31705
## [47,] 67.7094527 30.019905 74.160019 104.17992 104.59728
## [48,] 48.4991239 2.671690 56.596682 59.26837 59.63474
## [49,] 24.3928827 30.964743 98.078651 129.04339 129.53492
## [50,] 76.5459788 74.211966 57.681274 131.89324 132.21284cor(list)## y x1 x2 z1 z2
## y 1.0000000000 0.07758827 -0.07522322 -0.0005374415 -0.0008533985
## x1 0.0775882672 1.00000000 0.25313221 0.7815275201 0.7811302343
## x2 -0.0752232228 0.25313221 1.00000000 0.8013821453 0.8017562132
## z1 -0.0005374415 0.78152752 0.80138215 1.0000000000 0.9999955888
## z2 -0.0008533985 0.78113023 0.80175621 0.9999955888 1.0000000000#model catches exact multicollinearity easily
mcmodel1 <- lm(y ~ x1+x2+z1+z2)
summary(mcmodel1)##
## Call:
## lm(formula = y ~ x1 + x2 + z1 + z2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -48.767 -21.484 -2.126 19.999 53.375
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 51.22 10.53 4.865 1.38e-05 ***
## x1 15.81 28.05 0.563 0.576
## x2 15.65 28.09 0.557 0.580
## z1 NA NA NA NA
## z2 -15.72 28.06 -0.560 0.578
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 26.93 on 46 degrees of freedom
## Multiple R-squared: 0.0223, Adjusted R-squared: -0.04146
## F-statistic: 0.3498 on 3 and 46 DF, p-value: 0.7895#has a harder time catching near multicollinearity useful to use VIF or tolerance
mcmodel2 <- lm(y ~ x1+x2+z2)
summary(mcmodel2)##
## Call:
## lm(formula = y ~ x1 + x2 + z2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -48.767 -21.484 -2.126 19.999 53.375
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 51.22 10.53 4.865 1.38e-05 ***
## x1 15.81 28.05 0.563 0.576
## x2 15.65 28.09 0.557 0.580
## z2 -15.72 28.06 -0.560 0.578
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 26.93 on 46 degrees of freedom
## Multiple R-squared: 0.0223, Adjusted R-squared: -0.04146
## F-statistic: 0.3498 on 3 and 46 DF, p-value: 0.7895car::vif(mcmodel2)## x1 x2 z2
## 45305.13 49446.24 118711.40#Can drop the collinear term
mcmodel3 <- lm(y ~ x1+x2)
summary(mcmodel3)##
## Call:
## lm(formula = y ~ x1 + x2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -50.044 -19.925 -0.354 19.215 55.525
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 47.74591 8.44508 5.654 8.96e-07 ***
## x1 0.09332 0.13522 0.690 0.493
## x2 -0.08784 0.12964 -0.678 0.501
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 26.73 on 47 degrees of freedom
## Multiple R-squared: 0.01564, Adjusted R-squared: -0.02625
## F-statistic: 0.3733 on 2 and 47 DF, p-value: 0.6905car::vif(mcmodel3)## x1 x2
## 1.068463 1.068463#Or drop the other one
mcmodel4 <- lm(y ~ x1+z2)
summary(mcmodel4)##
## Call:
## lm(formula = y ~ x1 + z2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -50.050 -19.928 -0.363 19.217 55.519
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 47.77939 8.46296 5.646 9.21e-07 ***
## x1 0.18145 0.20950 0.866 0.391
## z2 -0.08807 0.12948 -0.680 0.500
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 26.73 on 47 degrees of freedom
## Multiple R-squared: 0.01571, Adjusted R-squared: -0.02618
## F-statistic: 0.3751 on 2 and 47 DF, p-value: 0.6893car::vif(mcmodel4)## x1 z2
## 2.565184 2.565184#normality check
#https://www.statmethods.net/stats/regression.html
nmodel <- summary(model2)
nmodel$residuals## 1 2 3 4 5 6
## 4.493754 55.950733 55.033809 -106.927120 -51.073841 37.578930
## 7 8 9 10 11 12
## 65.493754 120.893095 71.029094 -10.744020 -86.203964 -23.565690
## 13 14 15 16 17 18
## -75.687837 -84.819765 31.434310 -48.823026 43.515067 9.922899
## 19 20 21 22 23 24
## -64.562429 79.035614 25.236418 19.622717 -7.162278 43.016698
## 25 26 27 28 29 30
## -23.694358 -53.011137 -98.070580 68.556457 25.217205 -10.982143
## 31 32 33 34 35 36
## 11.029094 -44.321396 38.869976 102.559718 44.664107 -36.757062
## 37 38 39 40 41 42
## 17.874692 -107.851374 23.022573 72.302382 62.061348 -26.430205
## 43 44 45 46 47 48
## -54.710310 123.869976 -152.020273 -54.580187 10.826485 -122.311614
## 49 50 51 52 53 54
## 33.186756 -16.998095 -30.288845 -37.739008 20.415393 37.746024
## 55 56 57 58 59 60
## -74.819765 14.551742 129.285270 24.309199 -22.697618 4.248650layout(matrix(c(1,2,3,4),2,2)) # optional 4 graphs/page
plot(model2)
hist(nmodel$residuals)
#Runs different tests for normality run from the predicted values for rent
p3 <- predict(model2)
summary(p3)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 417.3 491.6 526.1 572.3 640.7 1045.8sd(p3)## [1] 125.7441#H0: from normal distribution p < 0.05 reject
shapiro.test(p3)##
## Shapiro-Wilk normality test
##
## data: p3
## W = 0.86842, p-value = 1.114e-05#compares if two samples are from same distribution so comparing to a normal distribution H0: from different distributions p < 0.05 reject
ks.test(p3, rnorm(60,572.3,125.7441))## Warning in ks.test(p3, rnorm(60, 572.3, 125.7441)): cannot compute exact p-value
## with ties##
## Two-sample Kolmogorov-Smirnov test
##
## data: p3 and rnorm(60, 572.3, 125.7441)
## D = 0.18333, p-value = 0.2656
## alternative hypothesis: two-sided
#Linearity Test
#https://bookdown.org/ccolonescu/RPoE4/further-inference-in-multiple-regression.html
#http://r-statistics.co/Statistical-Tests-in-R.html
library(lmtest)## Warning: package 'lmtest' was built under R version 4.0.2#Ramsey RESET test test whether higher order polynomials are necessary H0: no higher order polynomials are necssary
resettest(model3, power = 2:3, type = 'fitted')##
## RESET test
##
## data: model3
## RESET = 1.3333, df1 = 2, df2 = 52, p-value = 0.2725resettest(model3, power = 2:3, type = 'regressor')##
## RESET test
##
## data: model3
## RESET = 2.0443, df1 = 10, df2 = 44, p-value = 0.05111#Heteroscedasticity
#Put simply, heteroscedasticity (also spelled heteroskedasticity) refers to the circumstance in which the variability of a variable is unequal across the range of values of a second variable that predicts it.
#http://www.statsmakemecry.com/smmctheblog/confusing-stats-terms-explained-heteroscedasticity-heteroske.html
#Fisher Test can be used to tell if two samples have the same variance H0: ratio of variances is 1 aka they are the same p < 0.05 reject H0
var.test(sample(35,p3, replace = TRUE),sample(35,p3, replace = TRUE))##
## F test to compare two variances
##
## data: sample(35, p3, replace = TRUE) and sample(35, p3, replace = TRUE)
## F = 0.91901, num df = 513, denom df = 513, p-value = 0.3391
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.7728094 1.0928623
## sample estimates:
## ratio of variances
## 0.9190072#Independance
#Chi square tests if two caterogical variables are dependant on each other H0: variables are independant p < 0.05 reject H0
chisq.test(table(data1$Bedrooms,data1$Bathrooms))## Warning in chisq.test(table(data1$Bedrooms, data1$Bathrooms)): Chi-squared
## approximation may be incorrect##
## Pearson's Chi-squared test
##
## data: table(data1$Bedrooms, data1$Bathrooms)
## X-squared = 29.391, df = 3, p-value = 1.853e-06summary(table(data1$Bedrooms,data1$Bathrooms))## Number of cases in table: 60
## Number of factors: 2
## Test for independence of all factors:
## Chisq = 29.391, df = 3, p-value = 1.853e-06
## Chi-squared approximation may be incorrect#Correlation test between two continuous variables H0: correlation is 0 aka they are independant p < 0.05 reject H0
#All show some correlation
cor.test(data1$Rent, data1$Area)##
## Pearson's product-moment correlation
##
## data: data1$Rent and data1$Area
## t = 13.702, df = 58, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.7970327 0.9231055
## sample estimates:
## cor
## 0.8740597cor.test(data1$Bedrooms, data1$Area)##
## Pearson's product-moment correlation
##
## data: data1$Bedrooms and data1$Area
## t = 12.811, df = 58, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.7747757 0.9140118
## sample estimates:
## cor
## 0.8595894cor.test(data1$Bedrooms, data1$Bathrooms)##
## Pearson's product-moment correlation
##
## data: data1$Bedrooms and data1$Bathrooms
## t = 6.6217, df = 58, p-value = 1.262e-08
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.4826241 0.7800925
## sample estimates:
## cor
## 0.6561389cor.test(data1$Bedrooms, data1$Bathrooms)##
## Pearson's product-moment correlation
##
## data: data1$Bedrooms and data1$Bathrooms
## t = 6.6217, df = 58, p-value = 1.262e-08
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.4826241 0.7800925
## sample estimates:
## cor
## 0.6561389#Are uncorrelated
cor.test(data1$Area, data1$Distance)##
## Pearson's product-moment correlation
##
## data: data1$Area and data1$Distance
## t = -0.4492, df = 58, p-value = 0.655
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.3081967 0.1980052
## sample estimates:
## cor
## -0.05887988#Both of these tests use log parameters as well as lag and leads to determine if the variance changes or if the predictors are truly independant of each other#Other helpful Stat and R learning
#http://faculty.marshall.usc.edu/gareth-james/ISL/index.html
#https://web.stanford.edu/~hastie/ElemStatLearn/
#Essentials of Machine Learning https://www.analyticsvidhya.com/blog/2017/09/common-machine-learning-algorithms/
#CLT and Stat Basics https://www.analyticsvidhya.com/blog/2019/05/statistics-101-introduction-central-limit-theorem/Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I.
When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).
The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.
Erick Jones
PhD Candidate
Erick Jones is a Ph.D. candidate in Operations Research and Industrial Engineering who develops multi-systems optimization models to analyze how energy systems, water resources, supply chains, urban space, and transportation networks operate in concert to influence economic and environmental well-being.