Data 605 Final
Saayed Alam
5/16/2019
Problem 1
Using R, generate a random variable X that has 10,000 random uniform numbers from 1 to N, where N can be any number of your choosing greater than or equal to 6. Then generate a random variable Y that has 10,000 random normal numbers with a mean of \(\mu =\sigma = \frac{(N+1)}{2}\).
# load libraries
library(kableExtra)
library(psych)
library(dplyr)
library(ggplot2)
library(tidyr)
library(corrplot)
library(RColorBrewer)
library(matrixcalc)
library(MASS)
library(ggpubr)
library(gmodels)
library(mice)
library(e1071)
library(randomForest)
# generate random numbers
set.seed(123)
n <- 6
X <- runif(10000, 1, n)
Y <- rnorm(10000, (n+1)/2, (n+1)/2)
df <- data.frame(cbind(X, Y))
total <- nrow(df)Probability. Calculate as a minimum the below probabilities a through c. Assume the small letter “x” is estimated as the median of the X variable, and the small letter “y” is estimated as the 1st quartile of the Y variable. Interpret the meaning of all probabilities.
# "x" is estimated as the median of the X variable
x <- median(X)
# "y" is estimated as the 1st quartile of the Y variable
y <- summary(Y)[2]- \(P(X>x|X>y)\)
# P(X>y)
pX_gy <- nrow(subset(df, X > y))/total
# P(X>x & X>y)
pXgxXgy <- nrow(subset(df, X > x & Y > y))/total
# P(X>x | X>y)
p1 <- round(pXgxXgy / pX_gy, 4)
print(paste0("P(X>x | Y>y) = ", p1))## [1] "P(X>x | Y>y) = 0.3896"- \(P(X>x, Y>y)\)
# P(X>x & Y>y)
pXgxYgy <- nrow(subset(df, X > x & Y > y))/total
# P(X>x, Y>y)
print(paste0("P(X>x, Y>y) = ", pXgxYgy))## [1] "P(X>x, Y>y) = 0.3756"- \(P(X<x | X>y)\)
# P(X<x & X>y)
pXlxXgy <- nrow(subset(df, X < x & X > y))/total
# P(X<x | X>y)
p2 <- round(pXlxXgy / pX_gy, 4)
print(paste0("P(X<x | X>y) = ", p2))## [1] "P(X<x | X>y) = 0.4814"Investigate whether \(P(X>x \quad \& \quad Y>y)=P(X>x)P(Y>y)\) by building a table and evaluating the marginal and joint probabilities.
# marginal probablity of P(X>x)
pXgx <- nrow(subset(df, X > x))/total
# marginal probablity of P(Y>y)
pYgy <- nrow(subset(df, Y > y))/total
# joint probability of P(X>x & Y>y)
pXgxYgy <- nrow(subset(df, X > x & Y > y))/total
# product of marginal probabilities
product <- pXgx*pYgy
# joint probablity = product of marginal probablity
equal <- if(round(product, 2) == round(pXgxYgy, 2)){
print("Yes")
} else {
print("No")
}## [1] "Yes"# table
kable(cbind(pXgx, pYgy, product, equal, pXgxYgy), col.names = c("P(X>x)", "P(Y>y)", "P(X>x)P(Y>y)", "=", "P(X>x & Y>y)")) %>%
kable_styling("striped", full_width = F)| P(X>x) | P(Y>y) | P(X>x)P(Y>y) | = | P(X>x & Y>y) |
|---|---|---|---|---|
| 0.5 | 0.75 | 0.375 | Yes | 0.3756 |
Check to see if independence holds by using Fisher’s Exact Test and the Chi Square Test. What is the difference between the two? Which is most appropriate?
The difference between Fisher’s Exact Test and the Chi Square Test is Fisher’s Exact Test works well with small sample size. When Fisher’s Exact test is used with large sample size, result is an approximation. Since, we are working with large sample size of 10,000, it is more appropriate to use Chi Square test.
# contingency table
g_x <- subset(df, X > x)
g_y <- subset(df, Y > y)
le_x <- subset(df, X <= x)
le_y <- subset(df, Y <= y)
table <- matrix(c(nrow(g_x), nrow(g_y), nrow(le_x), nrow(le_y)), 2, 2,
dimnames = list(c("x", "y"),
c("X > x, Y > y", "X <= x, Y <= y")))
kable(table) %>%
kable_styling("striped", full_width = F)| X > x, Y > y | X <= x, Y <= y | |
|---|---|---|
| x | 5000 | 5000 |
| y | 7500 | 2500 |
# chi square test
chisq.test(table)##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: table
## X-squared = 1332.3, df = 1, p-value < 2.2e-16# fisher's exact test
fisher.test(table)##
## Fisher's Exact Test for Count Data
##
## data: table
## p-value < 2.2e-16
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.3137986 0.3540659
## sample estimates:
## odds ratio
## 0.3333333\(H_o\): X and Y are independent.
\(H_a\): X and Y are not independent.
Since the P-value is less than the significance level (0.05) for both test, we cannot accept the null hypothesis. Thus, we conclude that there is a relationship between X and Y.
Problem 2
Provide univariate descriptive statistics and appropriate plots for the training data set. Provide a scatterplot matrix for at least two of the independent variables and the dependent variable. Derive a correlation matrix for any three quantitative variables in the dataset. Test the hypotheses that the correlations between each pairwise set of variables is 0 and provide an 80% confidence interval. Discuss the meaning of your analysis. Would you be worried about familywise error? Why or why not?
# load data
hp_test <- read.csv("https://raw.githubusercontent.com/saayedalam/Data/master/house_prices_test.csv")
hp_train <- read.csv("https://raw.githubusercontent.com/saayedalam/Data/master/house_prices_train.csv")In the descriptive statistics below, we see our dataset has 1460 observations and 81 variables. We also see most of the variables are categorical and the rest have continuous values. Lastly, we see several of the variables have Na values as a factor.
# univariate descriptive statistics
glimpse(hp_train)## Observations: 1,460
## Variables: 81
## $ Id <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 1...
## $ MSSubClass <int> 60, 20, 60, 70, 60, 50, 20, 60, 50, 190, 20, 60,...
## $ MSZoning <fct> RL, RL, RL, RL, RL, RL, RL, RL, RM, RL, RL, RL, ...
## $ LotFrontage <int> 65, 80, 68, 60, 84, 85, 75, NA, 51, 50, 70, 85, ...
## $ LotArea <int> 8450, 9600, 11250, 9550, 14260, 14115, 10084, 10...
## $ Street <fct> Pave, Pave, Pave, Pave, Pave, Pave, Pave, Pave, ...
## $ Alley <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, ...
## $ LotShape <fct> Reg, Reg, IR1, IR1, IR1, IR1, Reg, IR1, Reg, Reg...
## $ LandContour <fct> Lvl, Lvl, Lvl, Lvl, Lvl, Lvl, Lvl, Lvl, Lvl, Lvl...
## $ Utilities <fct> AllPub, AllPub, AllPub, AllPub, AllPub, AllPub, ...
## $ LotConfig <fct> Inside, FR2, Inside, Corner, FR2, Inside, Inside...
## $ LandSlope <fct> Gtl, Gtl, Gtl, Gtl, Gtl, Gtl, Gtl, Gtl, Gtl, Gtl...
## $ Neighborhood <fct> CollgCr, Veenker, CollgCr, Crawfor, NoRidge, Mit...
## $ Condition1 <fct> Norm, Feedr, Norm, Norm, Norm, Norm, Norm, PosN,...
## $ Condition2 <fct> Norm, Norm, Norm, Norm, Norm, Norm, Norm, Norm, ...
## $ BldgType <fct> 1Fam, 1Fam, 1Fam, 1Fam, 1Fam, 1Fam, 1Fam, 1Fam, ...
## $ HouseStyle <fct> 2Story, 1Story, 2Story, 2Story, 2Story, 1.5Fin, ...
## $ OverallQual <int> 7, 6, 7, 7, 8, 5, 8, 7, 7, 5, 5, 9, 5, 7, 6, 7, ...
## $ OverallCond <int> 5, 8, 5, 5, 5, 5, 5, 6, 5, 6, 5, 5, 6, 5, 5, 8, ...
## $ YearBuilt <int> 2003, 1976, 2001, 1915, 2000, 1993, 2004, 1973, ...
## $ YearRemodAdd <int> 2003, 1976, 2002, 1970, 2000, 1995, 2005, 1973, ...
## $ RoofStyle <fct> Gable, Gable, Gable, Gable, Gable, Gable, Gable,...
## $ RoofMatl <fct> CompShg, CompShg, CompShg, CompShg, CompShg, Com...
## $ Exterior1st <fct> VinylSd, MetalSd, VinylSd, Wd Sdng, VinylSd, Vin...
## $ Exterior2nd <fct> VinylSd, MetalSd, VinylSd, Wd Shng, VinylSd, Vin...
## $ MasVnrType <fct> BrkFace, None, BrkFace, None, BrkFace, None, Sto...
## $ MasVnrArea <int> 196, 0, 162, 0, 350, 0, 186, 240, 0, 0, 0, 286, ...
## $ ExterQual <fct> Gd, TA, Gd, TA, Gd, TA, Gd, TA, TA, TA, TA, Ex, ...
## $ ExterCond <fct> TA, TA, TA, TA, TA, TA, TA, TA, TA, TA, TA, TA, ...
## $ Foundation <fct> PConc, CBlock, PConc, BrkTil, PConc, Wood, PConc...
## $ BsmtQual <fct> Gd, Gd, Gd, TA, Gd, Gd, Ex, Gd, TA, TA, TA, Ex, ...
## $ BsmtCond <fct> TA, TA, TA, Gd, TA, TA, TA, TA, TA, TA, TA, TA, ...
## $ BsmtExposure <fct> No, Gd, Mn, No, Av, No, Av, Mn, No, No, No, No, ...
## $ BsmtFinType1 <fct> GLQ, ALQ, GLQ, ALQ, GLQ, GLQ, GLQ, ALQ, Unf, GLQ...
## $ BsmtFinSF1 <int> 706, 978, 486, 216, 655, 732, 1369, 859, 0, 851,...
## $ BsmtFinType2 <fct> Unf, Unf, Unf, Unf, Unf, Unf, Unf, BLQ, Unf, Unf...
## $ BsmtFinSF2 <int> 0, 0, 0, 0, 0, 0, 0, 32, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ BsmtUnfSF <int> 150, 284, 434, 540, 490, 64, 317, 216, 952, 140,...
## $ TotalBsmtSF <int> 856, 1262, 920, 756, 1145, 796, 1686, 1107, 952,...
## $ Heating <fct> GasA, GasA, GasA, GasA, GasA, GasA, GasA, GasA, ...
## $ HeatingQC <fct> Ex, Ex, Ex, Gd, Ex, Ex, Ex, Ex, Gd, Ex, Ex, Ex, ...
## $ CentralAir <fct> Y, Y, Y, Y, Y, Y, Y, Y, Y, Y, Y, Y, Y, Y, Y, Y, ...
## $ Electrical <fct> SBrkr, SBrkr, SBrkr, SBrkr, SBrkr, SBrkr, SBrkr,...
## $ X1stFlrSF <int> 856, 1262, 920, 961, 1145, 796, 1694, 1107, 1022...
## $ X2ndFlrSF <int> 854, 0, 866, 756, 1053, 566, 0, 983, 752, 0, 0, ...
## $ LowQualFinSF <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ GrLivArea <int> 1710, 1262, 1786, 1717, 2198, 1362, 1694, 2090, ...
## $ BsmtFullBath <int> 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 0, ...
## $ BsmtHalfBath <int> 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ FullBath <int> 2, 2, 2, 1, 2, 1, 2, 2, 2, 1, 1, 3, 1, 2, 1, 1, ...
## $ HalfBath <int> 1, 0, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, ...
## $ BedroomAbvGr <int> 3, 3, 3, 3, 4, 1, 3, 3, 2, 2, 3, 4, 2, 3, 2, 2, ...
## $ KitchenAbvGr <int> 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 1, 1, 1, 1, 1, 1, ...
## $ KitchenQual <fct> Gd, TA, Gd, Gd, Gd, TA, Gd, TA, TA, TA, TA, Ex, ...
## $ TotRmsAbvGrd <int> 8, 6, 6, 7, 9, 5, 7, 7, 8, 5, 5, 11, 4, 7, 5, 5,...
## $ Functional <fct> Typ, Typ, Typ, Typ, Typ, Typ, Typ, Typ, Min1, Ty...
## $ Fireplaces <int> 0, 1, 1, 1, 1, 0, 1, 2, 2, 2, 0, 2, 0, 1, 1, 0, ...
## $ FireplaceQu <fct> NA, TA, TA, Gd, TA, NA, Gd, TA, TA, TA, NA, Gd, ...
## $ GarageType <fct> Attchd, Attchd, Attchd, Detchd, Attchd, Attchd, ...
## $ GarageYrBlt <int> 2003, 1976, 2001, 1998, 2000, 1993, 2004, 1973, ...
## $ GarageFinish <fct> RFn, RFn, RFn, Unf, RFn, Unf, RFn, RFn, Unf, RFn...
## $ GarageCars <int> 2, 2, 2, 3, 3, 2, 2, 2, 2, 1, 1, 3, 1, 3, 1, 2, ...
## $ GarageArea <int> 548, 460, 608, 642, 836, 480, 636, 484, 468, 205...
## $ GarageQual <fct> TA, TA, TA, TA, TA, TA, TA, TA, Fa, Gd, TA, TA, ...
## $ GarageCond <fct> TA, TA, TA, TA, TA, TA, TA, TA, TA, TA, TA, TA, ...
## $ PavedDrive <fct> Y, Y, Y, Y, Y, Y, Y, Y, Y, Y, Y, Y, Y, Y, Y, Y, ...
## $ WoodDeckSF <int> 0, 298, 0, 0, 192, 40, 255, 235, 90, 0, 0, 147, ...
## $ OpenPorchSF <int> 61, 0, 42, 35, 84, 30, 57, 204, 0, 4, 0, 21, 0, ...
## $ EnclosedPorch <int> 0, 0, 0, 272, 0, 0, 0, 228, 205, 0, 0, 0, 0, 0, ...
## $ X3SsnPorch <int> 0, 0, 0, 0, 0, 320, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...
## $ ScreenPorch <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 176, 0, 0, 0...
## $ PoolArea <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ PoolQC <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, ...
## $ Fence <fct> NA, NA, NA, NA, NA, MnPrv, NA, NA, NA, NA, NA, N...
## $ MiscFeature <fct> NA, NA, NA, NA, NA, Shed, NA, Shed, NA, NA, NA, ...
## $ MiscVal <int> 0, 0, 0, 0, 0, 700, 0, 350, 0, 0, 0, 0, 0, 0, 0,...
## $ MoSold <int> 2, 5, 9, 2, 12, 10, 8, 11, 4, 1, 2, 7, 9, 8, 5, ...
## $ YrSold <int> 2008, 2007, 2008, 2006, 2008, 2009, 2007, 2009, ...
## $ SaleType <fct> WD, WD, WD, WD, WD, WD, WD, WD, WD, WD, WD, New,...
## $ SaleCondition <fct> Normal, Normal, Normal, Abnorml, Normal, Normal,...
## $ SalePrice <int> 208500, 181500, 223500, 140000, 250000, 143000, ...summary(hp_train)## Id MSSubClass MSZoning LotFrontage
## Min. : 1.0 Min. : 20.0 C (all): 10 Min. : 21.00
## 1st Qu.: 365.8 1st Qu.: 20.0 FV : 65 1st Qu.: 59.00
## Median : 730.5 Median : 50.0 RH : 16 Median : 69.00
## Mean : 730.5 Mean : 56.9 RL :1151 Mean : 70.05
## 3rd Qu.:1095.2 3rd Qu.: 70.0 RM : 218 3rd Qu.: 80.00
## Max. :1460.0 Max. :190.0 Max. :313.00
## NA's :259
## LotArea Street Alley LotShape LandContour
## Min. : 1300 Grvl: 6 Grvl: 50 IR1:484 Bnk: 63
## 1st Qu.: 7554 Pave:1454 Pave: 41 IR2: 41 HLS: 50
## Median : 9478 NA's:1369 IR3: 10 Low: 36
## Mean : 10517 Reg:925 Lvl:1311
## 3rd Qu.: 11602
## Max. :215245
##
## Utilities LotConfig LandSlope Neighborhood Condition1
## AllPub:1459 Corner : 263 Gtl:1382 NAmes :225 Norm :1260
## NoSeWa: 1 CulDSac: 94 Mod: 65 CollgCr:150 Feedr : 81
## FR2 : 47 Sev: 13 OldTown:113 Artery : 48
## FR3 : 4 Edwards:100 RRAn : 26
## Inside :1052 Somerst: 86 PosN : 19
## Gilbert: 79 RRAe : 11
## (Other):707 (Other): 15
## Condition2 BldgType HouseStyle OverallQual
## Norm :1445 1Fam :1220 1Story :726 Min. : 1.000
## Feedr : 6 2fmCon: 31 2Story :445 1st Qu.: 5.000
## Artery : 2 Duplex: 52 1.5Fin :154 Median : 6.000
## PosN : 2 Twnhs : 43 SLvl : 65 Mean : 6.099
## RRNn : 2 TwnhsE: 114 SFoyer : 37 3rd Qu.: 7.000
## PosA : 1 1.5Unf : 14 Max. :10.000
## (Other): 2 (Other): 19
## OverallCond YearBuilt YearRemodAdd RoofStyle
## Min. :1.000 Min. :1872 Min. :1950 Flat : 13
## 1st Qu.:5.000 1st Qu.:1954 1st Qu.:1967 Gable :1141
## Median :5.000 Median :1973 Median :1994 Gambrel: 11
## Mean :5.575 Mean :1971 Mean :1985 Hip : 286
## 3rd Qu.:6.000 3rd Qu.:2000 3rd Qu.:2004 Mansard: 7
## Max. :9.000 Max. :2010 Max. :2010 Shed : 2
##
## RoofMatl Exterior1st Exterior2nd MasVnrType MasVnrArea
## CompShg:1434 VinylSd:515 VinylSd:504 BrkCmn : 15 Min. : 0.0
## Tar&Grv: 11 HdBoard:222 MetalSd:214 BrkFace:445 1st Qu.: 0.0
## WdShngl: 6 MetalSd:220 HdBoard:207 None :864 Median : 0.0
## WdShake: 5 Wd Sdng:206 Wd Sdng:197 Stone :128 Mean : 103.7
## ClyTile: 1 Plywood:108 Plywood:142 NA's : 8 3rd Qu.: 166.0
## Membran: 1 CemntBd: 61 CmentBd: 60 Max. :1600.0
## (Other): 2 (Other):128 (Other):136 NA's :8
## ExterQual ExterCond Foundation BsmtQual BsmtCond BsmtExposure
## Ex: 52 Ex: 3 BrkTil:146 Ex :121 Fa : 45 Av :221
## Fa: 14 Fa: 28 CBlock:634 Fa : 35 Gd : 65 Gd :134
## Gd:488 Gd: 146 PConc :647 Gd :618 Po : 2 Mn :114
## TA:906 Po: 1 Slab : 24 TA :649 TA :1311 No :953
## TA:1282 Stone : 6 NA's: 37 NA's: 37 NA's: 38
## Wood : 3
##
## BsmtFinType1 BsmtFinSF1 BsmtFinType2 BsmtFinSF2
## ALQ :220 Min. : 0.0 ALQ : 19 Min. : 0.00
## BLQ :148 1st Qu.: 0.0 BLQ : 33 1st Qu.: 0.00
## GLQ :418 Median : 383.5 GLQ : 14 Median : 0.00
## LwQ : 74 Mean : 443.6 LwQ : 46 Mean : 46.55
## Rec :133 3rd Qu.: 712.2 Rec : 54 3rd Qu.: 0.00
## Unf :430 Max. :5644.0 Unf :1256 Max. :1474.00
## NA's: 37 NA's: 38
## BsmtUnfSF TotalBsmtSF Heating HeatingQC CentralAir
## Min. : 0.0 Min. : 0.0 Floor: 1 Ex:741 N: 95
## 1st Qu.: 223.0 1st Qu.: 795.8 GasA :1428 Fa: 49 Y:1365
## Median : 477.5 Median : 991.5 GasW : 18 Gd:241
## Mean : 567.2 Mean :1057.4 Grav : 7 Po: 1
## 3rd Qu.: 808.0 3rd Qu.:1298.2 OthW : 2 TA:428
## Max. :2336.0 Max. :6110.0 Wall : 4
##
## Electrical X1stFlrSF X2ndFlrSF LowQualFinSF
## FuseA: 94 Min. : 334 Min. : 0 Min. : 0.000
## FuseF: 27 1st Qu.: 882 1st Qu.: 0 1st Qu.: 0.000
## FuseP: 3 Median :1087 Median : 0 Median : 0.000
## Mix : 1 Mean :1163 Mean : 347 Mean : 5.845
## SBrkr:1334 3rd Qu.:1391 3rd Qu.: 728 3rd Qu.: 0.000
## NA's : 1 Max. :4692 Max. :2065 Max. :572.000
##
## GrLivArea BsmtFullBath BsmtHalfBath FullBath
## Min. : 334 Min. :0.0000 Min. :0.00000 Min. :0.000
## 1st Qu.:1130 1st Qu.:0.0000 1st Qu.:0.00000 1st Qu.:1.000
## Median :1464 Median :0.0000 Median :0.00000 Median :2.000
## Mean :1515 Mean :0.4253 Mean :0.05753 Mean :1.565
## 3rd Qu.:1777 3rd Qu.:1.0000 3rd Qu.:0.00000 3rd Qu.:2.000
## Max. :5642 Max. :3.0000 Max. :2.00000 Max. :3.000
##
## HalfBath BedroomAbvGr KitchenAbvGr KitchenQual
## Min. :0.0000 Min. :0.000 Min. :0.000 Ex:100
## 1st Qu.:0.0000 1st Qu.:2.000 1st Qu.:1.000 Fa: 39
## Median :0.0000 Median :3.000 Median :1.000 Gd:586
## Mean :0.3829 Mean :2.866 Mean :1.047 TA:735
## 3rd Qu.:1.0000 3rd Qu.:3.000 3rd Qu.:1.000
## Max. :2.0000 Max. :8.000 Max. :3.000
##
## TotRmsAbvGrd Functional Fireplaces FireplaceQu GarageType
## Min. : 2.000 Maj1: 14 Min. :0.000 Ex : 24 2Types : 6
## 1st Qu.: 5.000 Maj2: 5 1st Qu.:0.000 Fa : 33 Attchd :870
## Median : 6.000 Min1: 31 Median :1.000 Gd :380 Basment: 19
## Mean : 6.518 Min2: 34 Mean :0.613 Po : 20 BuiltIn: 88
## 3rd Qu.: 7.000 Mod : 15 3rd Qu.:1.000 TA :313 CarPort: 9
## Max. :14.000 Sev : 1 Max. :3.000 NA's:690 Detchd :387
## Typ :1360 NA's : 81
## GarageYrBlt GarageFinish GarageCars GarageArea GarageQual
## Min. :1900 Fin :352 Min. :0.000 Min. : 0.0 Ex : 3
## 1st Qu.:1961 RFn :422 1st Qu.:1.000 1st Qu.: 334.5 Fa : 48
## Median :1980 Unf :605 Median :2.000 Median : 480.0 Gd : 14
## Mean :1979 NA's: 81 Mean :1.767 Mean : 473.0 Po : 3
## 3rd Qu.:2002 3rd Qu.:2.000 3rd Qu.: 576.0 TA :1311
## Max. :2010 Max. :4.000 Max. :1418.0 NA's: 81
## NA's :81
## GarageCond PavedDrive WoodDeckSF OpenPorchSF EnclosedPorch
## Ex : 2 N: 90 Min. : 0.00 Min. : 0.00 Min. : 0.00
## Fa : 35 P: 30 1st Qu.: 0.00 1st Qu.: 0.00 1st Qu.: 0.00
## Gd : 9 Y:1340 Median : 0.00 Median : 25.00 Median : 0.00
## Po : 7 Mean : 94.24 Mean : 46.66 Mean : 21.95
## TA :1326 3rd Qu.:168.00 3rd Qu.: 68.00 3rd Qu.: 0.00
## NA's: 81 Max. :857.00 Max. :547.00 Max. :552.00
##
## X3SsnPorch ScreenPorch PoolArea PoolQC
## Min. : 0.00 Min. : 0.00 Min. : 0.000 Ex : 2
## 1st Qu.: 0.00 1st Qu.: 0.00 1st Qu.: 0.000 Fa : 2
## Median : 0.00 Median : 0.00 Median : 0.000 Gd : 3
## Mean : 3.41 Mean : 15.06 Mean : 2.759 NA's:1453
## 3rd Qu.: 0.00 3rd Qu.: 0.00 3rd Qu.: 0.000
## Max. :508.00 Max. :480.00 Max. :738.000
##
## Fence MiscFeature MiscVal MoSold
## GdPrv: 59 Gar2: 2 Min. : 0.00 Min. : 1.000
## GdWo : 54 Othr: 2 1st Qu.: 0.00 1st Qu.: 5.000
## MnPrv: 157 Shed: 49 Median : 0.00 Median : 6.000
## MnWw : 11 TenC: 1 Mean : 43.49 Mean : 6.322
## NA's :1179 NA's:1406 3rd Qu.: 0.00 3rd Qu.: 8.000
## Max. :15500.00 Max. :12.000
##
## YrSold SaleType SaleCondition SalePrice
## Min. :2006 WD :1267 Abnorml: 101 Min. : 34900
## 1st Qu.:2007 New : 122 AdjLand: 4 1st Qu.:129975
## Median :2008 COD : 43 Alloca : 12 Median :163000
## Mean :2008 ConLD : 9 Family : 20 Mean :180921
## 3rd Qu.:2009 ConLI : 5 Normal :1198 3rd Qu.:214000
## Max. :2010 ConLw : 5 Partial: 125 Max. :755000
## (Other): 9In our pairs plot below using the pysch package, we see the histogram, scatterplot and correlation coefficient of all our variables. Since this a very a wide dataset, I have separated 10 variables per panel.
# plots for the training data set
hp_train %>%
dplyr::select(2:11) %>%
pairs.panels(method = "pearson", hist.col = "#f44542")
hp_train %>%
dplyr::select(12:21) %>%
pairs.panels(method = "pearson", hist.col = "#414df4")
hp_train %>%
dplyr::select(22:31) %>%
pairs.panels(method = "pearson", hist.col = "#c325fc")
hp_train %>%
dplyr::select(32:41) %>%
pairs.panels(method = "pearson", hist.col = "#ed0000")
hp_train %>%
dplyr::select(42:51) %>%
pairs.panels(method = "pearson", hist.col = "#00AFBB")
hp_train %>%
dplyr::select(52:61) %>%
pairs.panels(method = "pearson", hist.col = "#008233")
hp_train %>%
dplyr::select(62:71) %>%
pairs.panels(method = "pearson", hist.col = "#823800")
hp_train %>%
dplyr::select(72:81) %>%
pairs.panels(method = "pearson", hist.col = "#0004ff")
In our scatterplot matrices below, we picked two independent variables of our interest and the dependent variable. We selected LotArea, X1stFlrSF as our independent variables and our dependent variable is SalePrice. We can tell from the graph below, there is a weak to moderate correlation between these three variables. We can also tell all the variables are normally distributed with a right skewness.
#a scatterplot matrix of independent variables and the dependent variable
hp_train %>%
dplyr::select(c("SalePrice", "LotArea", "X1stFlrSF")) %>%
pairs.panels(method = "pearson", hist.col = "#c95656")
Preliminary test to check the test assumptions - Is the covariation linear? Yes, form the plot above, the relationship is linear. And the data looks normally distributed.
# correlation test between SalePrice and LotArea variables
pair1 <- cor.test(hp_train$SalePrice, hp_train$LotArea, method = "pearson", conf.level = 0.80)
pair2 <- cor.test(hp_train$SalePrice, hp_train$X1stFlrSF, method = "pearson", conf.level = 0.80)The p-value of the pair1 is 1.123139e-24, which is less than the significance level alpha = 0.05. We can conclude that SalePrice and LotArea are significantly correlated with a correlation coefficient of 0.2638434. There is sufficient evidence to conclude there is a significant linear relationship between SalePrice and LotArea because the correlation coefficient is significantly different from zero.
pair1##
## Pearson's product-moment correlation
##
## data: hp_train$SalePrice and hp_train$LotArea
## t = 10.445, df = 1458, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 80 percent confidence interval:
## 0.2323391 0.2947946
## sample estimates:
## cor
## 0.2638434The p-value of the pair2 is 5.394711e-147, which is less than the significance level alpha = 0.05. We can conclude that SalePrice and X1stFlrSF are significantly correlated with a correlation coefficient of 0.6058522. There is sufficient evidence to conclude there is a significant linear relationship between SalePrice and X1stFlrSF because the correlation coefficient is significantly different from zero.
pair2##
## Pearson's product-moment correlation
##
## data: hp_train$SalePrice and hp_train$X1stFlrSF
## t = 29.078, df = 1458, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 80 percent confidence interval:
## 0.5841687 0.6266715
## sample estimates:
## cor
## 0.6058522I would not be worried about familywise error because in both cases we rejected the null hypothesis. The familywise error is the probability of a coming to at least one false conclusion in a series of hypothesis tests. In both cases where the alternative hypothesis was affirmed and the p-values were extremely small.
# correlation plot
hp_corr <- hp_train %>%
dplyr::select(c("SalePrice", "LotArea", "X1stFlrSF")) %>%
cor()
corrplot.mixed(hp_corr,
lower = "number",
upper = "color",
lower.col = brewer.pal(n = 3, name = "Dark2"),
upper.col = brewer.pal(n = 3, name = "Dark2"))
Invert your correlation matrix from above. (This is known as the precision matrix and contains variance inflation factors on the diagonal.) Multiply the correlation matrix by the precision matrix, and then multiply the precision matrix by the correlation matrix. Conduct LU decomposition on the matrix.
# correlation matrix
hp_corr## SalePrice LotArea X1stFlrSF
## SalePrice 1.0000000 0.2638434 0.6058522
## LotArea 0.2638434 1.0000000 0.2994746
## X1stFlrSF 0.6058522 0.2994746 1.0000000# precision matrix
solve(hp_corr)## SalePrice LotArea X1stFlrSF
## SalePrice 1.5987636 -0.1447277 -0.9252722
## LotArea -0.1447277 1.1116223 -0.2452191
## X1stFlrSF -0.9252722 -0.2452191 1.6340150# multiply the correlation matrix by the precision matrix
hp_corr%*%solve(hp_corr)## SalePrice LotArea X1stFlrSF
## SalePrice 1 0 0
## LotArea 0 1 0
## X1stFlrSF 0 0 1# multiply the precision matrix by the correlation matrix
solve(hp_corr)%*%hp_corr## SalePrice LotArea X1stFlrSF
## SalePrice 1.000000e+00 0.000000e+00 0.000000e+00
## LotArea 2.775558e-17 1.000000e+00 2.775558e-17
## X1stFlrSF 0.000000e+00 -5.551115e-17 1.000000e+00# conduct LU decomposition on the matrix
luA <- lu.decomposition(hp_corr)
# test to show the product of LU is the original matrix
luA$L%*%luA$U## [,1] [,2] [,3]
## [1,] 1.0000000 0.2638434 0.6058522
## [2,] 0.2638434 1.0000000 0.2994746
## [3,] 0.6058522 0.2994746 1.0000000Many times, it makes sense to fit a closed form distribution to data. Select a variable in the Kaggle.com training dataset that is skewed to the right, shift it so that the minimum value is absolutely above zero if necessary. Then load the MASS package and run fitdistr to fit an exponential probability density function. Find the optimal value of \(\lambda\) for this distribution, and then take 1000 samples from this exponential distribution using this value (e.g., rexp(1000, \(\lambda\))). Plot a histogram and compare it with a histogram of your original variable. Using the exponential pdf, find the 5th and 95th percentiles using the cumulative distribution function (CDF). Also generate a 95% confidence interval from the empirical data, assuming normality. Finally, provide the empirical 5th percentile and 95th percentile of the data. Discuss.
We have selected the variable X1stFlrSF. The minimum value of our variable is 334; which is absolutely above zero.
# plot of the selected variable
ggplot(hp_train, aes(X1stFlrSF)) +
geom_histogram(aes(y = ..density..), colour="black", fill="#ed0000") +
geom_density(alpha=.5, fill="#f44542") +
geom_vline(xintercept = min(hp_train$X1stFlrSF)) +
labs(title = "First Floor in square feet") +
theme_minimal()
Then we load the MASS package and run fitdistr to fit an exponential probability density function. The optimal value of \(\lambda\) for this distribution is 0.0008601213
# fit exponential distribution
fit_dest <- fitdistr(hp_train$X1stFlrSF, "exponential")
rate <- fit_dest$estimate
rate## rate
## 0.0008601213In our histogram plots below, we take 1000 samples from the above exponential distribution using the value (e.g., rexp(1000, \(\lambda\))). Then, we plot a histogram and compare it with a histogram of your original variable. As expected, the exponential graph looks like an exponential distribution. However, our selected variable is normally distributed and skewed to the right.
# plot the exponential distribution
set.seed(123)
exp_lambda <- rexp(1000, rate)
exp_data <- as.data.frame(exp_lambda)
# transform the selected variable for plot
org_variable <- hp_train$X1stFlrSF
org_data <- as.data.frame(org_variable)
# histogram plot of exponential distribution
ed <- ggplot(exp_data, aes(exp_lambda)) +
geom_histogram(binwidth = function(x) 2 * IQR(x) / (length(x)^(1/3))) +
labs(x = "Exponential Distribution") +
theme_minimal()
# histogram plot of selected variable
ov <- ggplot(org_data, aes(org_variable)) +
geom_histogram(binwidth = function(x) 2 * IQR(x) / (length(x)^(1/3))) +
labs(x = "Selected Variable - 1st Floor in sq/ft") +
theme_minimal()
# comparing both plot
ggarrange(ed, ov, nrow = 1, ncol = 2)
Here, using the exponential pdf, we find the 5th and 95th percentiles are 61.129, 3627.486 respectively, using the cumulative distribution function (CDF). We also generate a 95% confidence interval 1142.78, 1182.47 from the empirical data, assuming normality. Finally, the empirical 5th percentile and 95th percentile of the data is 672.95, 1831.25, respectively. The selected variable is the X1stFlrSF from the empirical data. The empirical percentile mean that 95% of the 1460 house’s first floor are under 1,831 sq/ft. We also notice large disparity in the percentiles between the empirical data and the data fitted with exponential function. This suggests that the exponential function is not a good fit for this data set. The 95% confident interval built above was the confident interval of the mean. This means that if we sample the houses in that area repeatedly, each time with the same sample size (1460), then 95% of the time the mean of the sample will be between 1142.78, 1182.47. However, for a data this skewed, median is a much better descriptor than mean.
# the 5th and 95th percentiles of exponential pdf
exp_perc <- quantile(ecdf(exp_lambda), c(0.05, 0.95))
exp_perc## 5% 95%
## 61.12906 3627.48628# 95% confidence interval from the original data, assuming normality
ci(org_variable, confidence = 0.95)## Estimate CI lower CI upper Std. Error
## 1162.62671 1142.78038 1182.47304 10.11746# the 5th and 95th percentiles of original data
org_perc <- quantile(org_variable, c(0.05, 0.95))
org_perc## 5% 95%
## 672.95 1831.25Build some type of multiple regression model and submit your model to the competition board. Provide your complete model summary and results with analysis. Report your Kaggle.com user name and score.
First, we will preprocess our dataset before building a model. We will combine both train and test dataset. This will help us keep similar data transformation on both dataset. Later, we will separate them in to two different dataset to test our model. Second, we will split the full dataset based on their variable type - character and integer variables. Character variable’s Na values will replaced with Not Available as factor. All variable’s missing values will be imputed using the mice package. Lastly, we will combine both variable types and split them into train and test dataset.
# combine both dataset
full <- bind_rows(hp_train, hp_test)
# separate character variables
chr <- full[, sapply(full, is.character)]
# replace Na with a factor variable
chr[is.na(chr)] <- "Not Available"
fac <- chr %>%
lapply(as.factor) %>%
as.data.frame()
# separate integer variables
int <- full[, sapply(full, is.integer)]
# combine both variable type
full <- bind_cols(fac, int)
# impute the missing values
micemod <- full %>%
mice(method = "rf")##
## iter imp variable
## 1 1 LotFrontage MasVnrArea BsmtFinSF1 BsmtFinSF2 BsmtUnfSF TotalBsmtSF BsmtFullBath BsmtHalfBath GarageYrBlt GarageCars GarageArea SalePrice
## 1 2 LotFrontage MasVnrArea BsmtFinSF1 BsmtFinSF2 BsmtUnfSF TotalBsmtSF BsmtFullBath BsmtHalfBath GarageYrBlt GarageCars GarageArea SalePrice
## 1 3 LotFrontage MasVnrArea BsmtFinSF1 BsmtFinSF2 BsmtUnfSF TotalBsmtSF BsmtFullBath BsmtHalfBath GarageYrBlt GarageCars GarageArea SalePrice
## 1 4 LotFrontage MasVnrArea BsmtFinSF1 BsmtFinSF2 BsmtUnfSF TotalBsmtSF BsmtFullBath BsmtHalfBath GarageYrBlt GarageCars GarageArea SalePrice
## 1 5 LotFrontage MasVnrArea BsmtFinSF1 BsmtFinSF2 BsmtUnfSF TotalBsmtSF BsmtFullBath BsmtHalfBath GarageYrBlt GarageCars GarageArea SalePrice
## 2 1 LotFrontage MasVnrArea BsmtFinSF1 BsmtFinSF2 BsmtUnfSF TotalBsmtSF BsmtFullBath BsmtHalfBath GarageYrBlt GarageCars GarageArea SalePrice
## 2 2 LotFrontage MasVnrArea BsmtFinSF1 BsmtFinSF2 BsmtUnfSF TotalBsmtSF BsmtFullBath BsmtHalfBath GarageYrBlt GarageCars GarageArea SalePrice
## 2 3 LotFrontage MasVnrArea BsmtFinSF1 BsmtFinSF2 BsmtUnfSF TotalBsmtSF BsmtFullBath BsmtHalfBath GarageYrBlt GarageCars GarageArea SalePrice
## 2 4 LotFrontage MasVnrArea BsmtFinSF1 BsmtFinSF2 BsmtUnfSF TotalBsmtSF BsmtFullBath BsmtHalfBath GarageYrBlt GarageCars GarageArea SalePrice
## 2 5 LotFrontage MasVnrArea BsmtFinSF1 BsmtFinSF2 BsmtUnfSF TotalBsmtSF BsmtFullBath BsmtHalfBath GarageYrBlt GarageCars GarageArea SalePrice
## 3 1 LotFrontage MasVnrArea BsmtFinSF1 BsmtFinSF2 BsmtUnfSF TotalBsmtSF BsmtFullBath BsmtHalfBath GarageYrBlt GarageCars GarageArea SalePrice
## 3 2 LotFrontage MasVnrArea BsmtFinSF1 BsmtFinSF2 BsmtUnfSF TotalBsmtSF BsmtFullBath BsmtHalfBath GarageYrBlt GarageCars GarageArea SalePrice
## 3 3 LotFrontage MasVnrArea BsmtFinSF1 BsmtFinSF2 BsmtUnfSF TotalBsmtSF BsmtFullBath BsmtHalfBath GarageYrBlt GarageCars GarageArea SalePrice
## 3 4 LotFrontage MasVnrArea BsmtFinSF1 BsmtFinSF2 BsmtUnfSF TotalBsmtSF BsmtFullBath BsmtHalfBath GarageYrBlt GarageCars GarageArea SalePrice
## 3 5 LotFrontage MasVnrArea BsmtFinSF1 BsmtFinSF2 BsmtUnfSF TotalBsmtSF BsmtFullBath BsmtHalfBath GarageYrBlt GarageCars GarageArea SalePrice
## 4 1 LotFrontage MasVnrArea BsmtFinSF1 BsmtFinSF2 BsmtUnfSF TotalBsmtSF BsmtFullBath BsmtHalfBath GarageYrBlt GarageCars GarageArea SalePrice
## 4 2 LotFrontage MasVnrArea BsmtFinSF1 BsmtFinSF2 BsmtUnfSF TotalBsmtSF BsmtFullBath BsmtHalfBath GarageYrBlt GarageCars GarageArea SalePrice
## 4 3 LotFrontage MasVnrArea BsmtFinSF1 BsmtFinSF2 BsmtUnfSF TotalBsmtSF BsmtFullBath BsmtHalfBath GarageYrBlt GarageCars GarageArea SalePrice
## 4 4 LotFrontage MasVnrArea BsmtFinSF1 BsmtFinSF2 BsmtUnfSF TotalBsmtSF BsmtFullBath BsmtHalfBath GarageYrBlt GarageCars GarageArea SalePrice
## 4 5 LotFrontage MasVnrArea BsmtFinSF1 BsmtFinSF2 BsmtUnfSF TotalBsmtSF BsmtFullBath BsmtHalfBath GarageYrBlt GarageCars GarageArea SalePrice
## 5 1 LotFrontage MasVnrArea BsmtFinSF1 BsmtFinSF2 BsmtUnfSF TotalBsmtSF BsmtFullBath BsmtHalfBath GarageYrBlt GarageCars GarageArea SalePrice
## 5 2 LotFrontage MasVnrArea BsmtFinSF1 BsmtFinSF2 BsmtUnfSF TotalBsmtSF BsmtFullBath BsmtHalfBath GarageYrBlt GarageCars GarageArea SalePrice
## 5 3 LotFrontage MasVnrArea BsmtFinSF1 BsmtFinSF2 BsmtUnfSF TotalBsmtSF BsmtFullBath BsmtHalfBath GarageYrBlt GarageCars GarageArea SalePrice
## 5 4 LotFrontage MasVnrArea BsmtFinSF1 BsmtFinSF2 BsmtUnfSF TotalBsmtSF BsmtFullBath BsmtHalfBath GarageYrBlt GarageCars GarageArea SalePrice
## 5 5 LotFrontage MasVnrArea BsmtFinSF1 BsmtFinSF2 BsmtUnfSF TotalBsmtSF BsmtFullBath BsmtHalfBath GarageYrBlt GarageCars GarageArea SalePricefull <- complete(micemod)
# split train and test dataset
train <- full[1:length(hp_train$SalePrice),]
test <- full[(length(hp_train$SalePrice)+1):nrow(full),]Finally, we use two different algorithms to see which yields the better result. We use SVM from e1071 package. And randomforest package. Random forest technique gave us the best result and kaggle score with username is below.
# svm model
svm_model <- svm(SalePrice ~ ., data = train, cost = 3)
svm_pred <- predict(svm_model, newdata = test)
# random forest model
rf_model <- randomForest(SalePrice ~ ., data = train)
rf_pred <- predict(rf_model, newdata = test)
# create submission file
submission <- as.data.frame(cbind(hp_test$Id, svm_pred))
sub1 <- as.data.frame(cbind(hp_test$Id, rf_pred))
colnames(submission) <- c("Id", "SalePrice")
colnames(sub1) <- c("Id", "SalePrice")
write.csv(submission, file = "Kaggle Competition 1", quote = FALSE, row.names = FALSE)
write.csv(sub1, file = "Kaggle Competition 1.1", quote = FALSE, row.names = FALSE)
summary(rf_pred)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 65739 130816 159103 179445 211175 489855# kaggle score
info <- c("saayedalam", 0.15250)
names(info) <- c("Username", "Public Score")
kable(info, col.names = "Kaggle") %>%
kable_styling(full_width = F)| Kaggle | |
|---|---|
| Username | saayedalam |
| Public Score | 0.1525 |
Link to my Kaggle account and Youtube video: 
