Housing price data from Ames, Iowa. 
The full dataset is available at https://ww2.amstat.org/publications/jse/v19n3/decock/AmesHousing.txt
A description of the full dataset is at https://ww2.amstat.org/publications/jse/v19n3/decock/DataDocumentation.txt 
We will be working with a slightly smaller version of the full dataset. 
Set the right working directory (you need to change the directory to the local directory in your machine containing the data)

```{r }
dataDir<-"../../finalDataSets"
```

Read the dataset:

```{r }
dd = read.csv(file.path(dataDir, "Ames_Short.csv"), header = T)
```

The size of the dataset (number of rows and columns) can be obtained by

```{r }
dim(dd)
```

Variable names:

```{r }
names(dd)
```



```{r }
head(dd)
```



```{r }
summary(dd)
```


Dataset can also be found here \url{https://archive.ics.uci.edu/ml/datasets/Bike+Sharing+Dataset} 
Read the data

```{r }
bike <- read.csv(file.path(dataDir, "BikeSharingDataset.csv"))
```



```{r }
dim(bike)
```



```{r }
names(bike)
```



```{r }
head(bike)
```



```{r }
summary(bike)
```

https://collegescorecard.ed.gov/assets/CollegeScorecardDataDictionary-10-26-2016.xlsx
https://collegescorecard.ed.gov/assets/FullDataDocumentation.pdf

```{r }
scorecard = read.csv(file.path(dataDir, "college_short.csv"))
```



```{r }
dim(scorecard)
```



```{r }
names(scorecard)
```



```{r }
head(scorecard)
```



```{r }
summary(scorecard)
```


Explanation of the dataset can be found at http://lib.stat.cmu.edu/datasets/bodyfat

```{r }
body = read.csv(file.path(dataDir, "bodyfat_short.csv"), header = T)
```



```{r }
dim(body)
```



```{r }
names(body)
```



```{r }
head(body)
```



```{r }
summary(body)
```


This is a dataset from the book "Introductory Econometrics: A Modern Approach" by Wooldridge. 

```{r }
load(file.path(dataDir, "campus.Rdata"))
```



```{r }
campus = data
```



```{r }
campus.desc = desc
```

This gives data for 97 colleges and universities for the year 1992. 

```{r }
#campus
```



```{r }
summary(campus)
```



```{r }
plot(crime ~ police, data = campus)
```



```{r }
m1 = lm(crime ~ police, data = campus)
summary(m1)
plot(crime ~ police, data = campus)
abline(m1)
```



```{r }
m2 = lm(lcrime ~ lpolice, data = campus)
summary(m2)
plot(lcrime ~ lpolice, data = campus)
abline(m2)
```



```{r }
y = campus$lcrime
x = campus$lpolice
nrep = 100
bslope = rep(-999, nrep)
bint = rep(-999, nrep)
for(i in 1:nrep)
    {
        
        sampled = sample(1:length(y), size = length(y), replace = T)
        bslope[i] = coef(lm(y[sampled] ~ x[sampled]))[2]
        bint[i] = coef(lm(y[sampled] ~ x[sampled]))[1]
        }
#quartz()
#par(mfrow = c(1, 2))
#hist(bslope, breaks = 500)
#abline(v = coef(m2)[2], col = "red")
#hist(bint, breaks = 500)
#abline(v = coef(m2)[1], col = "red")
#par(mfrow = c(1, 1))
```

The bootstrap standard errors are given by: 

```{r }
c(sd(bint), sd(bslope))
```

These seem slightly larger than the standard errors outputted by the lm() function:

```{r }
summary(m2)
```

Explanation of the dataset can be found at http://lib.stat.cmu.edu/datasets/bodyfat

```{r }
body <- read.csv(file.path(dataDir, "bodyfat_short.csv"), header = T)
```



```{r }
dim(body)
```



```{r }
names(body)
```



```{r }
head(body)
```



```{r }
summary(body)
```

Notice the unusual values in the bodyfat, age, weight and height variables.  

```{r }
body[body$HEIGHT < 30,]
```



```{r }
par(mfrow = c(3, 3))
for (i in 1:8)
    {
        hist(body[,i], xlab = "", main = names(body)[i], breaks = 500)
        }
par(mfrow = c(1, 1))
```

Again notice the presence of outliers. 

```{r }
body[body$HIP > 140,]
```

Subject 39 seems to have large values in many variables. 

```{r }
pairs(body)
```



```{r }
ou1 = which(body$HEIGHT < 30)
ou2 = which(body$WEIGHT > 300)
ou3 = which(body$HIP > 120)
ou = c(ou1, ou2, ou3)
pairs(body[-ou,])
```



```{r }
par(mfrow = c(3, 3))
for(i in 2:8)
{
    plot(body[,i], body[,1], xlab = names(body)[i], ylab = "BODYFAT")
    }
par(mfrow = c(1, 1))
```



```{r }
par(mfrow = c(3, 3))
for(i in 2:8)
{
    plot(body[-ou,i], body[-ou,1], xlab = names(body)[i], ylab = "BODYFAT")
    }
par(mfrow = c(1, 1))
```


Multiple Regression

```{r }
ft = lm(BODYFAT ~ AGE + WEIGHT + HEIGHT + CHEST + ABDOMEN + HIP + THIGH, data = body)
summary(ft)
```


Sum of squares of the regression equation given by R:

```{r }
eqn = -37.48 + 0.012*body$AGE - 0.139*body$WEIGHT - 0.102*body$HEIGHT - 0.0008*body$CHEST + 0.968*body$ABDOMEN - 0.183*body$HIP + 0.286*body$THIGH 
sum.sq = sum((body$BODYFAT - eqn)^2)
sum.sq
```


Sum of squares of my arbitrary equation:

```{r }
my.eqn = -30 + 0.05*body$AGE - 0.013*body$WEIGHT - 0.2*body$HEIGHT - 0.8*body$CHEST + 0.0009*body$ABDOMEN - 0.03*body$HIP + 0.6*body$THIGH 
my.sum.sq = sum((body$BODYFAT - my.eqn)^2)
my.sum.sq
```


Correlation between chest circumeference and abdomen circumference

```{r }
cor(body$CHEST, body$ABDOMEN)
```


Another regression for the bodyfat dataset (we dropped the variables Abdomen and Hip)

```{r }
ft1 = lm(BODYFAT ~ AGE + WEIGHT + HEIGHT + CHEST + THIGH, data = body) 
summary(ft1)
```


Regression for bodyfat percentage using only age and chest circumference

```{r }
ft2 = lm(BODYFAT ~ AGE + CHEST, data = body) 
summary(ft2)
```


Plotting the 3d scatterplot along with the fitted regression PLANE

```{r }
library(scatterplot3d)
sp = scatterplot3d(body$AGE, body$CHEST, body$BODYFAT)
sp$plane3d(ft2, lty.box = "solid", draw_lines = TRUE, draw_polygon = TRUE, col = "red")
```


Prediction using the regression equation

```{r }
bf.pred = -37.48 + 0.01202*30 - 0.1392*180 - 0.1028*70 - 0.0008312*90 + 0.9685*86 - 0.1834*97 + 0.2857*60
bf.pred
```


First observation (row) of the bodyfat dataset

```{r }
obs1 = body[1,]
obs1
```


Fitted value for the first obsevation

```{r }
-37.48 + 0.01202*body[1,"AGE"] - 0.1392*body[1,"WEIGHT"] - 0.1028*body[1,"HEIGHT"] - 0.0008312*body[1, "CHEST"] + 0.9685*body[1,"ABDOMEN"] - 0.1834*body[1, "HIP"] + 0.2857*body[1, "THIGH"]  
```


Fitted values in R

```{r }
ft = lm(BODYFAT ~ AGE + WEIGHT + HEIGHT + CHEST + ABDOMEN + HIP + THIGH, data = body)
names(ft)
head(ft$fitted.values)
```


Plotting fitted values against response

```{r }
plot(ft$fitted.values, body$BODYFAT, xlab = "Fitted Values", ylab = "Bodyfat Percentage")
```


Coefficient of Determination (or Multiple R^2)

```{r }
ft = lm(BODYFAT ~ AGE + WEIGHT + HEIGHT + CHEST + ABDOMEN + HIP + THIGH, data = body)
Rsquared = (cor(body$BODYFAT, ft$fitted.values))^2
Rsquared
```


R^2 in lm summary

```{r }
ft = lm(BODYFAT ~ AGE + WEIGHT + HEIGHT + CHEST + ABDOMEN + HIP + THIGH, data = body)
summary(ft)
```


Residual for the first observation

```{r }
ft = lm(BODYFAT ~ AGE + WEIGHT + HEIGHT + CHEST + ABDOMEN + HIP + THIGH, data = body)
res.1 = body[1,"BODYFAT"] - ft$fitted.values[1]
res.1
```


Residuals via lm

```{r }
ft = lm(BODYFAT ~ AGE + WEIGHT + HEIGHT + CHEST + ABDOMEN + HIP + THIGH, data = body)
names(ft)
head(ft$residuals)
```


Residuals plotted against the fitted values

```{r }
ft = lm(BODYFAT ~ AGE + WEIGHT + HEIGHT + CHEST + ABDOMEN + HIP + THIGH, data = body)
plot(ft$fitted.values, ft$residuals, xlab = "Fitted Values", ylab = "Residuals")
```


Residuals plotted against the explanatory variables

```{r }
par(mfrow = c(3, 3))
for(i in 2:8)
{
    plot(body[,i], ft$residuals, xlab = names(body)[i], ylab = "Residuals")
    }
par(mfrow = c(1, 1))
# Note here the presence of the outliers. 
```

Residuals are uncorrelated with every explanatory variable

```{r }
ft = lm(BODYFAT ~ AGE + WEIGHT + HEIGHT + CHEST + ABDOMEN + HIP + THIGH, data = body)
cor(ft$residuals, body$AGE)
cor(ft$residuals, body$WEIGHT)
cor(ft$residuals, body$HEIGHT)
cor(ft$residuals, body$CHEST)
cor(ft$residuals, body$ABDOMEN)
cor(ft$residuals, body$HIP)
cor(ft$residuals, body$THIGH)
```


Residuals have mean zero. 

```{r }
ft = lm(BODYFAT ~ AGE + WEIGHT + HEIGHT + CHEST + ABDOMEN + HIP + THIGH, data = body)
mean(ft$residuals)
```


Fitting a linear model to residuals:

```{r }
ft = lm(BODYFAT ~ AGE + WEIGHT + HEIGHT + CHEST + ABDOMEN + HIP + THIGH, data = body)
m.res = lm(ft$residuals ~ body$AGE + body$WEIGHT + body$HEIGHT + body$CHEST + body$ABDOMEN + body$HIP + body$THIGH)
summary(m.res)
```


Residual Sum of Squares

```{r }
ft = lm(BODYFAT ~ AGE + WEIGHT + HEIGHT + CHEST + ABDOMEN + HIP + THIGH, data = body)
rss.ft = sum((ft$residuals)^2)
rss.ft
```


Residual Sum of Squares and R^2

```{r }
ft = lm(BODYFAT ~ AGE + WEIGHT + HEIGHT + CHEST + ABDOMEN + HIP + THIGH, data = body)
rss.ft = sum((ft$residuals)^2)
rss.ft
tss = sum(((body$BODYFAT) - mean(body$BODYFAT))^2)
1 - (rss.ft/tss)
summary(ft)
```


RSS after removing the variable Abdomen

```{r }
ft.1 = lm(BODYFAT ~ AGE + WEIGHT + HEIGHT + CHEST + HIP + THIGH, data = body)
rss.ft1 = sum((ft.1$residuals)^2)
rss.ft1
```

Compare this to the RSS of the full regression. 

```{r }
rss.ft
```


RSS after removing the variable CHEST

```{r }
ft.2 = lm(BODYFAT ~ AGE + WEIGHT + HEIGHT + ABDOMEN + HIP + THIGH, data = body)
rss.ft2 = sum((ft.2$residuals)^2)
rss.ft2
```

Compare this to the RSS of the full regression. 

```{r }
rss.ft
```


R2 after removing the variable Abdomen

```{r }
ft.1 = lm(BODYFAT ~ AGE + WEIGHT + HEIGHT + CHEST + HIP + THIGH, data = body)
R2.ft1 = (cor(body$BODYFAT, ft.1$fitted.values))^2
R2.ft1
```

Compare this to the RSS of the full regression. 

```{r }
R2.ft  = (cor(body$BODYFAT, ft$fitted.values))^2
R2.ft
```


R2 after removing the variable chest

```{r }
ft.2 = lm(BODYFAT ~ AGE + WEIGHT + HEIGHT + ABDOMEN + HIP + THIGH, data = body)
R2.ft2 = (cor(body$BODYFAT, ft.2$fitted.values))^2
R2.ft2
```

Compare this to the RSS of the full regression. 

```{r }
R2.ft
```


Residual Degrees of Freedom

```{r }
ft = lm(BODYFAT ~ AGE + WEIGHT + HEIGHT + CHEST + ABDOMEN + HIP + THIGH, data = body)
n = nrow(body)
p = 7
rs.df = n-p-1
rs.df
```


Residual Standard Error

```{r }
ft = lm(BODYFAT ~ AGE + WEIGHT + HEIGHT + CHEST + ABDOMEN + HIP + THIGH, data = body)
rss = sum((ft$residuals)^2)
rse = sqrt(rss/rs.df)
rse
```


Summary (look for residual degrees of freedom and residual sum of squares)

```{r }
summary(ft)
```


https://collegescorecard.ed.gov/assets/CollegeScorecardDataDictionary-10-26-2016.xlsx
https://collegescorecard.ed.gov/assets/FullDataDocumentation.pdf
Set the right working directory

```{r }
scorecard <- read.csv(file.path(dataDir, "college_short.csv"))
```

Number of observations and Variables

```{r }
dim(scorecard)
```

Variable Names

```{r }
names(scorecard)
```


Fitting a regression equation in the usual way.  

```{r }
req.bad = lm(RET_FT4~ TUITIONFEE_OUT + TYPE, data = scorecard)
summary(req.bad)
```

How can the coefficients be interpreted here?
Is type being treated as a numeric variable by R?  

```{r }
is.numeric(scorecard$TYPE)
```


Regression equation:  

```{r }
req = lm(RET_FT4~ TUITIONFEE_OUT + as.factor(TYPE), data = scorecard)
summary(req)
```


Plotting the three regression lines in one plot:  

```{r }
plot(RET_FT4 ~ TUITIONFEE_OUT, data = scorecard, xlab = "Tuition Fee (out of state)", ylab = "Retention Rate", type = "n") 
w1 = (scorecard$TYPE == 1)
points(RET_FT4[w1] ~ TUITIONFEE_OUT[w1], data = scorecard, col = "blue")
abline(c(req$coefficients[1], req$coefficients[2]), col = "blue")
w2 = (scorecard$TYPE == 2)
points(RET_FT4[w2] ~ TUITIONFEE_OUT[w2], data = scorecard, col = "red")
abline(c(req$coefficients[1] + req$coefficients[3], req$coefficients[2]), col = "red")
w3 = (scorecard$TYPE == 3)
points(RET_FT4[w3] ~ TUITIONFEE_OUT[w3], data = scorecard, col = "green")
abline(c(req$coefficients[1] + req$coefficients[4], req$coefficients[2]), col = "green")
```


Regression equation with interaction:  

```{r }
req.1 = lm(RET_FT4~ TUITIONFEE_OUT + as.factor(TYPE) + TUITIONFEE_OUT : as.factor(TYPE), data = scorecard) 
summary(req.1)
```


Regression equation with interaction (another equivalent way):  

```{r }
req.2 = lm(RET_FT4~ TUITIONFEE_OUT * as.factor(TYPE), data = scorecard) 
summary(req.2)
```


Plotting the three regression lines in one plot:  

```{r }
plot(RET_FT4 ~ TUITIONFEE_OUT, data = scorecard, xlab = "Tuition Fee (out of state)", ylab = "Retention Rate", type = "n") 
w1 = (scorecard$TYPE == 1)
points(RET_FT4[w1] ~ TUITIONFEE_OUT[w1], data = scorecard, col = "blue")
abline(c(req.1$coefficients[1], req.1$coefficients[2]), col = "blue")
w2 = (scorecard$TYPE == 2)
points(RET_FT4[w2] ~ TUITIONFEE_OUT[w2], data = scorecard, col = "red")
abline(c(req.1$coefficients[1] + req.1$coefficients[3], req.1$coefficients[2] + req.1$coefficients[5]), col = "red")
w3 = (scorecard$TYPE == 3)
points(RET_FT4[w3] ~ TUITIONFEE_OUT[w3], data = scorecard, col = "green")
abline(c(req.1$coefficients[1] + req.1$coefficients[4], req.1$coefficients[2] + req.1$coefficients[6]), col = "green") 
```


Dataset can also be found here \url{https://archive.ics.uci.edu/ml/datasets/Bike+Sharing+Dataset} 
Set the right working directory
Read the data

```{r }
bike = read.csv(file.path(dataDir, "BikeSharingDataset.csv"))
```



```{r }
dim(bike)
```



```{r }
names(bike)
```



```{r }
head(bike)
```


Summary Statistics

```{r }
summary(bike$atemp)
summary(as.factor(bike$workingday))
summary(as.factor(bike$weathersit))
```



```{r }
#Regression Equation
md1 = lm(casual ~ atemp + as.factor(workingday) +
as.factor(weathersit), data = bike)
summary(md1)
```



```{r }
#Regression Equation
md2 = lm(casual ~ atemp + as.factor(workingday) +
as.factor(weathersit) + as.factor(workingday) : as.factor(weathersit),
data = bike) 
summary(md2)
```



```{r }
#Regression Equation
md3 = lm(casual ~ atemp + as.factor(workingday) +
as.factor(workingday) : atemp + as.factor(weathersit), data = bike)  
summary(md3)
```


This is a dataset from the book "Introductory Econometrics: A Modern Approach" by Wooldridge. 

```{r }
load(file.path(dataDir, "wage1.Rdata"))
wages = data
wages.desc = desc
```

This gives data on wage (hourly earning) and various other variables; data collected from 1976.  


```{r }
head(wages)
dim(wages)
wages.desc
```


This is a dataset from the book "Introductory Econometrics: A Modern Approach" by Wooldridge. 

```{r }
load(file.path(dataDir, "campus.Rdata"))
campus = data
campus.desc = desc
```

This gives data for 97 colleges and universities for the year 1992. 

```{r }
summary(campus)
```


This is a dataset from the book "Introductory Econometrics: A Modern Approach" by Wooldridge. 

```{r }
load(file.path(dataDir, "twoyear.Rdata"))
ty = data
ty.desc = desc
head(ty)
```

Data is collected for a sample of working people with a high school degree. One objective of analyzing this data is to comare the returns to education at junior colleges an four-year colleges (universities). 
This is a dataset from the book "Introductory Econometrics: A Modern Approach" by Wooldridge. 

```{r }
load(file.path(dataDir, "mlb1.Rdata"))
bb = data
bb.desc = desc
bb.desc
head(bb)
```


This is a dataset from the book "Introductory Econometrics: A Modern Approach" by Wooldridge. 

```{r }
load(file.path(dataDir, "wage1.Rdata"))
wages = data
wages.desc = desc
```

This gives data on wage (hourly earning) and various other variables; data collected from 1976.  

The model M

```{r }
M = lm(lwage ~ educ + exper + tenure, data = wages)
summary(M)
```


The model m

```{r }
m = lm(lwage ~ educ + tenure, data = wages)
summary(m)
```


This is a dataset from the book "Introductory Econometrics: A Modern Approach" by Wooldridge. 

```{r }
load(file.path(dataDir, "campus.Rdata"))
campus = data
campus.desc = desc
```

This gives data for 97 colleges and universities for the year 1992. 

The model M

```{r }
M = lm(lcrime ~ lenroll + lpolice, data = campus)
summary(M)
```


The model m

```{r }
m = lm(lcrime ~ offset(1*lenroll) + lpolice, data = campus)
summary(m)
```


This is a dataset from the book "Introductory Econometrics: A Modern Approach" by Wooldridge. 

```{r }
load(file.path(dataDir, "twoyear.Rdata"))
ty = data
ty.desc = desc
```

Data is collected for a sample of working people with a high school degree. One objective of analyzing this data is to comare the returns to education at junior colleges an four-year colleges (universities). 
The model M

```{r }
M = lm(lwage ~ jc + univ + exper, data = ty)
summary(M)
```


The model m

```{r }
m = lm(lwage ~ I(jc + univ) + exper, data = ty)
summary(m)
```


This is a dataset from the book "Introductory Econometrics: A Modern Approach" by Wooldridge. 

```{r }
load(file.path(dataDir, "mlb1.Rdata"))
bb = data
bb.desc = desc
```


The model M

```{r }
M = lm(lsalary ~ years + gamesyr + bavg + hrunsyr + rbisyr, data = bb)
summary(M)
```


The model m

```{r }
m = lm(lsalary ~ years + gamesyr, data = bb)
summary(m)
```


Plotting the density of the $F$-distribution
Choose first the degrees of freedom:

```{r }
d1 = 5
d2 = 30
xx = seq(0, 10, 0.01)
plot(xx, df(xx, d1, d2), type = "l", xlab = "x", ylab = "F-density")
```


This is a dataset from the book "Introductory Econometrics: A Modern Approach" by Wooldridge. 

```{r }
load(file.path(dataDir, "wage1.Rdata"))
wages = data
wages.desc = desc
```

This gives data on wage (hourly earning) and various other variables; data collected from 1976.  

The model M

```{r }
M = lm(lwage ~ educ + exper + tenure, data = wages)
summary(M)
```


The model m

```{r }
m = lm(lwage ~ educ + tenure, data = wages)
summary(m)
```


Calculating the F-statistic.

```{r }
rss.M = sum((M$residuals)^2)
rss.M
p = 3
rss.m = sum((m$residuals)^2)
rss.m
q = 2
n = nrow(wages)
F = ((rss.m - rss.M)/(p-q))/(rss.M/(n-p-1))
F
```


Plotting the density of the $F$-distribution

```{r }
d1 = p-q
d2 = n-p-1
xx = seq(0, 10, 0.01)
plot(xx, df(xx, d1, d2), type = "l", xlab = "x", ylab = "F-density")
abline(v = F)
```


Computing the p-value

```{r }
1 - pf(F, d1, d2)
```


p-value for testing H_0: beta_2 = 0 from summary:

```{r }
summary(M)
```


Testing via the anova function:

```{r }
anova(m, M)
```


The model M

```{r }
M = lm(lcrime ~ lenroll + lpolice, data = campus)
```


The model m

```{r }
m = lm(lcrime ~ offset(1*lenroll) + lpolice, data = campus)
```


Testing m against M

```{r }
anova(m, M)
```


The model M

```{r }
M = lm(lwage ~ jc + univ + exper, data = ty)
```


The submodel m

```{r }
m = lm(lwage ~ I(jc + univ) + exper, data = ty)
```


Test

```{r }
anova(m, M)
```




```{r }
M = lm(lsalary ~ years + gamesyr + bavg + hrunsyr + rbisyr, data = bb)
```




```{r }
m = lm(lsalary ~ years + gamesyr, data = bb)
```


Test

```{r }
anova(m, M)
```


Illustrating the F statistic given by R in summary()

```{r }
M = lm(lsalary ~ years + gamesyr + bavg + hrunsyr + rbisyr, data = bb)
summary(M)
```


Obtaining the F-statistic given by R summary via the anova function

```{r }
m = lm(lsalary ~ 1, data = bb)
summary(m)
anova(m, M)
```

Back to the bodyfat datset 

```{r }
body = read.csv(file.path(dataDir, "bodyfat_short.csv"), header = T)
#Let us fit a linear regression equation to bodyfat based on Age, Weight, Height, Chest, Abdomen, Hip and Thigh:
md = lm(BODYFAT ~ AGE + WEIGHT + HEIGHT + CHEST + ABDOMEN + HIP + THIGH, data = body)
summary(md)
```

Confidence interval

```{r }
n = nrow(body)
p = 7
t.cut.off = qt(0.975, (n-p-1))
t.cut.off
```

The 95% confidence interval for the true regression coefficient corresponding to ABDOMEN is: 

```{r }
lci = 0.9685 - t.cut.off*(0.08531)
uci = 0.9685 + t.cut.off*(0.08531)
c(lci, uci)
```


Confidence intervals for all the parameters:

```{r }
confint(md)
```


Prediction

```{r }
x0 = c(1, 30, 180, 70, 95, 90, 100, 60)
pred.bodyfat = sum(x0*md$coefficients)
pred.bodyfat
```


Intervals using the predict function

```{r }
x0 = data.frame(AGE = 30, WEIGHT = 180, HEIGHT = 70, CHEST = 95, ABDOMEN = 90, HIP = 100, THIGH = 60)
predict(md, x0, interval = "confidence")
predict(md, x0, interval = "prediction")
```


The car seat position dataset

```{r }
library(faraway)
data(seatpos)
```

help(seatpos)

```{r }
names(seatpos)
pairs(seatpos)
g = lm(hipcenter ~ ., seatpos)
summary(g)
```


Stepwise Procedures
We illustrate the backward method - at each stage we remove the
predictor with the largest p-value over a critical value (say, chosen to be 0.05):

```{r }
g <- update(g, . ~ . - Ht)
summary(g)
g <- update(g, . ~ . - Weight)
summary(g)
g <- update(g, . ~ . - Seated)
summary(g)
g <- update(g, . ~ . - Arm)
summary(g)
g <- update(g, . ~ . - Thigh)
summary(g)
g <- update(g, . ~ . - Age)
summary(g)
g <- update(g, . ~ . - Leg)
summary(g)
```

Forward Selection for the car seat position data

```{r }
for( i in 1:8)
  {
    g1 <- lm(hipcenter ~ ., seatpos[,c(i,9)])
    print((summary(g1))$coef)
}
```

The first selected predictor: Ht 

```{r }
for( i in c(1:3,5:8))
  {
     g2 <- lm(hipcenter ~ ., seatpos[,c(i,4,9)])
     print((summary(g2))$coef)
}
```

The second selected predictor: Leg
Continue this. 
Regsubsets on the car seat position data:

```{r }
library(leaps)
b = regsubsets(hipcenter ~ ., seatpos)
rs = summary(b)
rs$which
```


Leave One Out Cross Validation Score for model using all explanatory variables.

```{r }
n <- nrow(seatpos)
pred.y <- rep(NA,n) 
for (i in 1:nrow(seatpos))
  {
   g <- lm(hipcenter~., seatpos, subset=(1:n)[-i])
   pred.y[i] <- predict(g, seatpos[i,-9])
  }
cv.err.full <-  sum((seatpos[,9]-pred.y)^2)
cv.err.full
```



```{r }
for (i in 1:nrow(seatpos))
  {
   g <- lm(hipcenter ~ Age+Ht+Leg, seatpos, subset=(1:n)[-i])
   pred.y[i] <- predict(g, seatpos[i,c("Age","Ht","Leg")])
  }
cv.err.Age.Ht.Leg <-  sum((seatpos[,9]-pred.y)^2)
cv.err.Age.Ht.Leg
```

Back to the bodyfat datset 

```{r }
body = read.csv(file.path(dataDir, "bodyfat_short.csv"), header = T)
#Let us fit a linear regression equation to bodyfat based on Age, Weight, Height, Chest, Abdomen, Hip and Thigh:
md = lm(BODYFAT ~ AGE + WEIGHT + HEIGHT + CHEST + ABDOMEN + HIP + THIGH, data = body)
summary(md)
```

Regression diagnostics:

```{r }
par(mfrow = c(2, 2))
plot(md)
par(mfrow = c(1, 1))
```



```{r }
body[c(39, 42, 36),]
apply(body, 2, mean)
```

Non-linearity

```{r }
n = 200
xx = 3 + 4*abs(rnorm(n))
yy = -2 + 0.5*xx^(1.85) + rnorm(n)
m1 = lm(yy~xx)
par(mfrow = c(2, 2))
plot(m1)
par(mfrow = c(1, 1))
plot(yy~xx)
```


Including square of x as another explanatory variable

```{r }
m2 = lm(yy ~ xx + I(xx^2))
summary(m2)
par(mfrow = c(2, 2))
plot(m2)
par(mfrow = c(1, 1))
```


An example involving heteroscedasticity

```{r }
n = 200
xx = 3 + 4*abs(rnorm(n))
yy = -2 + 5*xx + (xx^(1.5))*rnorm(n)
m1 = lm(yy ~ xx)
par(mfrow = c(2, 2))
plot(m1)
par(mfrow = c(1, 1))
plot(yy~xx)
```


