Assignment_06
Mohamed Elashri
April 4, 2020
Question 1
first lets input the data
temp <- c(273,283,293,303,313,323,333,343,353,363,373)
vap_pres <- c(4.6,9.2,17.5,31.8,55.3,92.5,149.4,233.7,355.1,525.8,760.0)
data <- data.frame(temp,vap_pres)
rm(temp,vap_pres)a
we plot scatter plot
plot(x=data$temp, y=data$vap_pres , xlab="Tempreature", ylab="Vapor Preasure", main= "Vapor Preasure vs Tempreature ")
As we can see, this is curved graph so linear model will not fit the data and will not be adequate
b
we can fit a linear model and see the parameters.
linear <- lm(data$vap_pres~data$temp)
summary(linear)##
## Call:
## lm(formula = data$vap_pres ~ data$temp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -120.63 -92.10 -37.66 64.33 222.55
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1956.258 363.807 -5.377 0.000446 ***
## data$temp 6.686 1.121 5.964 0.000212 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 117.6 on 9 degrees of freedom
## Multiple R-squared: 0.7981, Adjusted R-squared: 0.7756
## F-statistic: 35.57 on 1 and 9 DF, p-value: 0.0002117so our model will be
y = 6.686x - 1956.258 B_0 = -1956.258 B_1 = 6.686
c
we can plot the residual as follows. also with other scal and leverage plots.
par (mfrow=c(2,2))
plot(linear)
from residuals vs fitted values plot, there is high variation in the residuals and they are far from zero line. this indicate that our model is not explaining the variation of preasure.
d
res = rstandard(linear)
qqnorm(res,
ylab="Preasure",
xlab="Tempreature",
main="Preasure explained vs Tempreature")
qqline(res)
As the model data follows our imagined normal distribution closely. we don’t think that there is problem with normality assumption. but there are two points that can be considered as outliers.
e
we can do log transformation
data$log_vap_pres <- log(data$vap_pres)
data$inv_tem <- -1/data$temp
non_linear = lm(log_vap_pres~inv_tem, data=data)
summary(non_linear)##
## Call:
## lm(formula = log_vap_pres ~ inv_tem, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.030999 -0.013107 0.004863 0.016735 0.021260
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.061e+01 6.325e-02 325.8 <2e-16 ***
## inv_tem 5.201e+03 2.014e+01 258.3 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02067 on 9 degrees of freedom
## Multiple R-squared: 0.9999, Adjusted R-squared: 0.9999
## F-statistic: 6.672e+04 on 1 and 9 DF, p-value: < 2.2e-16this model have R2 of 99.99% which indicate that this model is a very good model in explaining variation of preasure due to tempreature.
lets plot the residual
par (mfrow=c(2,2))
plot(non_linear)
residual now is closer to the zero line, this indicate that very small residuals and hence our model is good.
Question 2 (6.1)
first let us import our data.
data2 <- read.table("solar.csv", header= TRUE, sep=",", stringsAsFactors=FALSE, dec=".")then we proceed with head of our data to see variable names
head(data2)## y x1 x2 x3 x4 x5
## 1 271.8 783.35 33.53 40.55 16.66 13.20
## 2 264.0 748.45 36.50 36.19 16.46 14.11
## 3 238.8 684.45 34.66 37.31 17.66 15.68
## 4 230.7 827.80 33.13 32.52 17.50 10.53
## 5 251.6 860.45 35.75 33.71 16.40 11.00
## 6 257.9 875.15 34.46 34.14 16.28 11.31we fit the model now.
model2 <- lm(y ~ x1 + x2 + x3 + x4 + x5, data = data2)now let us now do some plotting
plot(model2)
We have now 4 graphs that gives us our analysis results but let us discuss it in more details.
** from residual vs fitted values plot, we can see that it indicate that there is no linear relations in the data.
** Scale function graph shows that there is a pattern with non horizontal line so we can confirm that there is no constant variance.
** from Residual vs Leverage graph, we can see that observations 1 and 4 are our high leverage points, and also they are with observation 22 the most influential points according to their cooks distance.
Question 3 (6.2)
we repeat everything from problem 2 and then discuss the results at the end
first let us import our data.
data3 <- read.table("property.csv", header= TRUE, sep=",", stringsAsFactors=FALSE, dec=".")then we proceed with head of our data to see variable names
head(data3)## y x1 x2 x3 x4 x5 x6 x7 x8 x9
## 1 29.5 5.0208 1 3.531 1.500 2 7 4 62 0
## 2 27.9 4.5429 1 2.275 1.175 1 6 3 40 0
## 3 25.9 4.5573 1 4.050 1.232 1 6 3 54 0
## 4 29.9 5.0597 1 4.455 1.121 1 6 3 42 0
## 5 29.9 3.8910 1 4.455 0.988 1 6 3 56 0
## 6 30.9 5.8980 1 5.850 1.240 1 7 3 51 1now procced with the creating a model
model3 <- lm(y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9, data = data3)now with plotting
plot(model3)
We have now 4 graphs that gives us our analysis results but let us discuss it in more details.
** from residual vs fitted values plot, we can see that that it suggests that there is linear relation.
** Scale function graph shows that there is a pattern with increasing horizontal line so we cannot reject constant variance assumption.
** from Residual vs Leverage graph, we can see that observations 16 and 14 are our high leverage points, and also they are with observation 7 the most influential points according to their cooks distance.
Question 4 (10.1)
As usual, we need to import our data first.
data4 <- read.table("football.csv", header= TRUE, sep=",", stringsAsFactors=FALSE, dec=".")lets explore data head
head(data4)## y x1 x2 x3 x4 x5 x6 x7 x8 x9
## 1 10 2113 1985 38.9 64.7 4 868 59.7 2205 1917
## 2 11 2003 2855 38.8 61.3 3 615 55.0 2096 1575
## 3 11 2957 1737 40.1 60.0 14 914 65.6 1847 2175
## 4 13 2285 2905 41.6 45.3 -4 957 61.4 1903 2476
## 5 10 2971 1666 39.2 53.8 15 836 66.1 1457 1866
## 6 11 2309 2927 39.7 74.1 8 786 61.0 1848 2339now let us construct our MLR model
model4 <- lm(y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9, data = data4)and using function to create forward, stepwise and backward directions
a
we use forward election algorithm as follows
Forward <- step(model4, direction=c("forward"))## Start: AIC=41.48
## y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9summary(Forward)##
## Call:
## lm(formula = y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9,
## data = data4)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.0408 -0.6802 -0.1131 0.9835 2.9785
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.292e+00 1.281e+01 -0.569 0.576312
## x1 8.124e-04 2.006e-03 0.405 0.690329
## x2 3.631e-03 8.410e-04 4.318 0.000414 ***
## x3 1.222e-01 2.590e-01 0.472 0.642750
## x4 3.189e-02 4.160e-02 0.767 0.453289
## x5 1.511e-05 4.684e-02 0.000 0.999746
## x6 1.590e-03 3.248e-03 0.490 0.630338
## x7 1.544e-01 1.521e-01 1.015 0.323547
## x8 -3.895e-03 2.052e-03 -1.898 0.073793 .
## x9 -1.791e-03 1.417e-03 -1.264 0.222490
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.83 on 18 degrees of freedom
## Multiple R-squared: 0.8156, Adjusted R-squared: 0.7234
## F-statistic: 8.846 on 9 and 18 DF, p-value: 5.303e-05b
then we use a backward elimination algorithm to select a subset regression model.
Backward <- step(model4, direction=c("backward"))## Start: AIC=41.48
## y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9
##
## Df Sum of Sq RSS AIC
## - x5 1 0.000 60.293 39.476
## - x1 1 0.549 60.842 39.730
## - x3 1 0.746 61.039 39.821
## - x6 1 0.803 61.096 39.847
## - x4 1 1.968 62.261 40.376
## - x7 1 3.451 63.744 41.035
## <none> 60.293 41.476
## - x9 1 5.348 65.642 41.856
## - x8 1 12.072 72.365 44.587
## - x2 1 62.448 122.741 59.380
##
## Step: AIC=39.48
## y ~ x1 + x2 + x3 + x4 + x6 + x7 + x8 + x9
##
## Df Sum of Sq RSS AIC
## - x1 1 0.553 60.846 37.732
## - x3 1 0.750 61.043 37.822
## - x6 1 0.818 61.111 37.854
## - x4 1 2.053 62.346 38.414
## - x7 1 3.859 64.152 39.213
## <none> 60.293 39.476
## - x9 1 5.351 65.644 39.857
## - x8 1 12.086 72.379 42.592
## - x2 1 66.979 127.272 58.395
##
## Step: AIC=37.73
## y ~ x2 + x3 + x4 + x6 + x7 + x8 + x9
##
## Df Sum of Sq RSS AIC
## - x6 1 0.690 61.536 36.048
## - x3 1 1.715 62.561 36.510
## - x4 1 3.051 63.897 37.102
## <none> 60.846 37.732
## - x9 1 4.852 65.698 37.880
## - x7 1 8.961 69.807 39.579
## - x8 1 16.599 77.445 42.486
## - x2 1 67.010 127.856 56.524
##
## Step: AIC=36.05
## y ~ x2 + x3 + x4 + x7 + x8 + x9
##
## Df Sum of Sq RSS AIC
## - x3 1 1.726 63.262 34.822
## - x4 1 2.767 64.303 35.279
## <none> 61.536 36.048
## - x9 1 4.831 66.367 36.164
## - x7 1 9.390 70.926 38.024
## - x8 1 18.314 79.851 41.343
## - x2 1 66.447 127.984 54.552
##
## Step: AIC=34.82
## y ~ x2 + x4 + x7 + x8 + x9
##
## Df Sum of Sq RSS AIC
## - x4 1 1.743 65.004 33.583
## <none> 63.262 34.822
## - x9 1 5.629 68.891 35.209
## - x8 1 17.701 80.962 39.730
## - x7 1 18.583 81.845 40.033
## - x2 1 75.598 138.860 54.835
##
## Step: AIC=33.58
## y ~ x2 + x7 + x8 + x9
##
## Df Sum of Sq RSS AIC
## <none> 65.004 33.583
## - x9 1 4.866 69.870 33.604
## - x7 1 16.908 81.913 38.057
## - x8 1 23.299 88.303 40.160
## - x2 1 82.892 147.897 54.601summary(Backward)##
## Call:
## lm(formula = y ~ x2 + x7 + x8 + x9, data = data4)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.3519 -0.5612 -0.0856 0.6972 3.2802
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.8217034 7.7847061 -0.234 0.81705
## x2 0.0038186 0.0007051 5.416 1.67e-05 ***
## x7 0.2168941 0.0886759 2.446 0.02252 *
## x8 -0.0040149 0.0013983 -2.871 0.00863 **
## x9 -0.0016349 0.0012460 -1.312 0.20244
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.681 on 23 degrees of freedom
## Multiple R-squared: 0.8012, Adjusted R-squared: 0.7666
## F-statistic: 23.17 on 4 and 23 DF, p-value: 8.735e-08c
Finally we use stepwise regression to select a subset regression model.
Stepwise <- step(model4, direction=c("both"))## Start: AIC=41.48
## y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9
##
## Df Sum of Sq RSS AIC
## - x5 1 0.000 60.293 39.476
## - x1 1 0.549 60.842 39.730
## - x3 1 0.746 61.039 39.821
## - x6 1 0.803 61.096 39.847
## - x4 1 1.968 62.261 40.376
## - x7 1 3.451 63.744 41.035
## <none> 60.293 41.476
## - x9 1 5.348 65.642 41.856
## - x8 1 12.072 72.365 44.587
## - x2 1 62.448 122.741 59.380
##
## Step: AIC=39.48
## y ~ x1 + x2 + x3 + x4 + x6 + x7 + x8 + x9
##
## Df Sum of Sq RSS AIC
## - x1 1 0.553 60.846 37.732
## - x3 1 0.750 61.043 37.822
## - x6 1 0.818 61.111 37.854
## - x4 1 2.053 62.346 38.414
## - x7 1 3.859 64.152 39.213
## <none> 60.293 39.476
## - x9 1 5.351 65.644 39.857
## + x5 1 0.000 60.293 41.476
## - x8 1 12.086 72.379 42.592
## - x2 1 66.979 127.272 58.395
##
## Step: AIC=37.73
## y ~ x2 + x3 + x4 + x6 + x7 + x8 + x9
##
## Df Sum of Sq RSS AIC
## - x6 1 0.690 61.536 36.048
## - x3 1 1.715 62.561 36.510
## - x4 1 3.051 63.897 37.102
## <none> 60.846 37.732
## - x9 1 4.852 65.698 37.880
## + x1 1 0.553 60.293 39.476
## - x7 1 8.961 69.807 39.579
## + x5 1 0.004 60.842 39.730
## - x8 1 16.599 77.445 42.486
## - x2 1 67.010 127.856 56.524
##
## Step: AIC=36.05
## y ~ x2 + x3 + x4 + x7 + x8 + x9
##
## Df Sum of Sq RSS AIC
## - x3 1 1.726 63.262 34.822
## - x4 1 2.767 64.303 35.279
## <none> 61.536 36.048
## - x9 1 4.831 66.367 36.164
## + x6 1 0.690 60.846 37.732
## + x1 1 0.425 61.111 37.854
## - x7 1 9.390 70.926 38.024
## + x5 1 0.027 61.509 38.036
## - x8 1 18.314 79.851 41.343
## - x2 1 66.447 127.984 54.552
##
## Step: AIC=34.82
## y ~ x2 + x4 + x7 + x8 + x9
##
## Df Sum of Sq RSS AIC
## - x4 1 1.743 65.004 33.583
## <none> 63.262 34.822
## - x9 1 5.629 68.891 35.209
## + x3 1 1.726 61.536 36.048
## + x1 1 1.321 61.941 36.231
## + x6 1 0.700 62.561 36.510
## + x5 1 0.101 63.160 36.777
## - x8 1 17.701 80.962 39.730
## - x7 1 18.583 81.845 40.033
## - x2 1 75.598 138.860 54.835
##
## Step: AIC=33.58
## y ~ x2 + x7 + x8 + x9
##
## Df Sum of Sq RSS AIC
## <none> 65.004 33.583
## - x9 1 4.866 69.870 33.604
## + x1 1 1.865 63.140 34.768
## + x4 1 1.743 63.262 34.822
## + x3 1 0.701 64.303 35.279
## + x6 1 0.451 64.554 35.388
## + x5 1 0.327 64.678 35.442
## - x7 1 16.908 81.913 38.057
## - x8 1 23.299 88.303 40.160
## - x2 1 82.892 147.897 54.601summary(Stepwise)##
## Call:
## lm(formula = y ~ x2 + x7 + x8 + x9, data = data4)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.3519 -0.5612 -0.0856 0.6972 3.2802
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.8217034 7.7847061 -0.234 0.81705
## x2 0.0038186 0.0007051 5.416 1.67e-05 ***
## x7 0.2168941 0.0886759 2.446 0.02252 *
## x8 -0.0040149 0.0013983 -2.871 0.00863 **
## x9 -0.0016349 0.0012460 -1.312 0.20244
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.681 on 23 degrees of freedom
## Multiple R-squared: 0.8012, Adjusted R-squared: 0.7666
## F-statistic: 23.17 on 4 and 23 DF, p-value: 8.735e-08d
for our three models we can see their value for Adjusted R-Squared.
** for forward, all the 9 variables get selected and Adjusted R-Square = 72.34% ** for backward, only the regressors x2,x7,x8 and x9 are selected and Adjusted R-Square = 76.66% ** for Stepwise, only the regressors x2,x7,x8 and x9 are selected and Adjusted R-Square = 76.66%
With backward and stepwise have the same and higher than forward Adjusted R-Square, we should use them as a method of selection.
Question 5 (10.2)
We are doing the same but now our model will contain only the assigned regressors and comapre it with the full model. I’m a litte bit confused with what all-possible-regressions procedures means. but what I understand is that it means that we should be considering the subset of all the assigned 6 individual regressors , but do not need to consider the interaction of them. so we will have 2^6 = 64 models.
Now let us define our model of 6 individual regressors.
using regsubsets from leaps package we can get the answer to what is the best subset based on our method (backward or forward).
library(leaps)
regsubsets <-regsubsets(y ~ x1 + x2 + x4 + x7 + x8 + x9,
data = data4,
nbest = 1, # 1 best model for each number of predictors
nvmax = 6, # NULL for no limit on number of variables
force.in = NULL, force.out = NULL,
method = c("backward","forward"))
regsubsets## Subset selection object
## Call: regsubsets.formula(y ~ x1 + x2 + x4 + x7 + x8 + x9, data = data4,
## nbest = 1, nvmax = 6, force.in = NULL, force.out = NULL,
## method = c("backward", "forward"))
## 6 Variables (and intercept)
## Forced in Forced out
## x1 FALSE FALSE
## x2 FALSE FALSE
## x4 FALSE FALSE
## x7 FALSE FALSE
## x8 FALSE FALSE
## x9 FALSE FALSE
## 1 subsets of each size up to 6
## Selection Algorithm: backwardnow lets get the summary of subsets.
summary(regsubsets)## Subset selection object
## Call: regsubsets.formula(y ~ x1 + x2 + x4 + x7 + x8 + x9, data = data4,
## nbest = 1, nvmax = 6, force.in = NULL, force.out = NULL,
## method = c("backward", "forward"))
## 6 Variables (and intercept)
## Forced in Forced out
## x1 FALSE FALSE
## x2 FALSE FALSE
## x4 FALSE FALSE
## x7 FALSE FALSE
## x8 FALSE FALSE
## x9 FALSE FALSE
## 1 subsets of each size up to 6
## Selection Algorithm: backward
## x1 x2 x4 x7 x8 x9
## 1 ( 1 ) " " " " " " " " "*" " "
## 2 ( 1 ) " " "*" " " " " "*" " "
## 3 ( 1 ) " " "*" " " "*" "*" " "
## 4 ( 1 ) " " "*" " " "*" "*" "*"
## 5 ( 1 ) "*" "*" " " "*" "*" "*"
## 6 ( 1 ) "*" "*" "*" "*" "*" "*"The function summary() reports the best set of variables for each model size. From the output above, an asterisk specifies that a given variable is included in the corresponding model.
As we can see, the best two variables model will contain x2 and x8. the best three variables model will contain x2, x7 and x8.
Now we want to know what is the best model from our 64 sub models we have. this will be done using Adjusted R^2, C_p and BIC values.
we get the best of these using following minimizing code.
res.sum <- summary(regsubsets)
data.frame(
Adj.R2 = which.max(res.sum$adjr2),
CP = which.min(res.sum$cp),
BIC = which.min(res.sum$bic)
)## Adj.R2 CP BIC
## 1 4 3 3adjusted R2 tells us that the best model is the one with 4 predictor variables. However, using the BIC and Cp criteria, we should go for the model with 3 variables.
Now lets plot RSS, adjusted R2, C_p and BIC vs number of variables we have
# Set up a 2x2 grid so we can look at 4 plots at once
par(mfrow = c(2,2))
plot(res.sum$rss, xlab = "Number of Variables", ylab = "RSS", type = "l")
plot(res.sum$adjr2, xlab = "Number of Variables", ylab = "Adjusted RSq", type = "l")
# We will now plot a red dot to indicate the model with the largest adjusted R^2 statistic.
# The which.max() function can be used to identify the location of the maximum point of a vector
adj_r2_max = which.max(res.sum$adjr2) # 11
# The points() command works like the plot() command, except that it puts points
# on a plot that has already been created instead of creating a new plot
points(adj_r2_max, res.sum$adjr2[adj_r2_max], col ="red", cex = 2, pch = 20)
# We'll do the same for C_p and BIC, this time looking for the models with the SMALLEST statistic
plot(res.sum$cp, xlab = "Number of Variables", ylab = "Cp", type = "l")
cp_min = which.min(res.sum$cp) # 10
points(cp_min, res.sum$cp[cp_min], col = "red", cex = 2, pch = 20)
plot(res.sum$bic, xlab = "Number of Variables", ylab = "BIC", type = "l")
bic_min = which.min(res.sum$bic) # 6
points(bic_min, res.sum$bic[bic_min], col = "red", cex = 2, pch = 20)
plot(regsubsets, scale = "adjr2")
Adjusted R2 is maximized by the model that contains x2, x7 and x8
coef(regsubsets, 3)## (Intercept) x2 x7 x8
## -1.808372059 0.003598070 0.193960210 -0.004815494plot(regsubsets, scale = "Cp")
Cp confirms this observation
coef(regsubsets, 3)## (Intercept) x2 x7 x8
## -1.808372059 0.003598070 0.193960210 -0.004815494So we can conclude that the best model of subsets will be the one with x2,x7 and x8 regressors.