01-ch1.Rmd

```{r setup, include=FALSE}
require(here)
options(scipen=99)
options(digits = 4)
```


# Linear Regression with One Predictor Variable {#ch1}


## Relations between Variables

### Example 1

```{r, fig.align='center', fig.height=5,fig.width=7,message=FALSE,warning=FALSE}

require(here)

fig1.2 <- read.csv(here('/data/fig1.2.csv/'),header=FALSE)

require(ggplot2)

ggplot(fig1.2,aes(x=V1,y=V2))+
  geom_point(size=4,col='gray42')+
  geom_smooth(method=lm,se=FALSE,lty=2,col='black')+
  xlim(55,100)+
  ylim(55,100)+
  xlab('Midyear Evaluation')+
  ylab('Year-End Evaluation')+
  theme_bw()
```


### Example 2

```{r, fig.align='center', fig.height=4,fig.width=7,message=FALSE,warning=FALSE}

fig1.3 <- read.csv(here('/data/fig1.3.csv/'),header=FALSE)

ggplot(fig1.3,aes(x=V1,y=V2))+
  geom_point(size=4,col='gray42')+
  geom_smooth(method=lm,formula = y ~ x + I(x^2), se=FALSE,lty=2,col='black')+
  xlim(7.5,25)+
  ylim(1,30)+
  xlab('Age')+
  ylab('Stereoid Level')+
  theme_bw()
```

## Regression Model an Their Uses

## Simple Linear Regression Model with Distributions of Error Terms Unspecified

### Example (page 10)

$$Y_i = 9.5 + 2.1X_i + \epsilon_i$$
```{r, fig.align='center', fig.height=7,fig.width=7,message=FALSE,warning=FALSE}

b0 = 9.5
b1 = 2.1
        
x = 0:70
        
mean = b0 + b1*x
    
err.sd = 5
        
  plot(x,mean,type="l",
       ylim=c(40,120),
       xlim=c(10,50),
       cex=1,
       pch=19,
       xlab="Number of Bids Prepared",
       ylab="Hours",
       xaxt='n')
        
  axis(side=1,at=c(0,25,45))
  abline(b0,b1)
        
  dens = dnorm(seq(-3,3,.01),0,1)

  for(i in c(26,46)){
            
    x. = x[i] - 5*dens
    y. = mean[i]+seq(-3,3,.01)*err.sd
    points(x.,y.,type="l",lty=2)
    abline(v=x[i],lty=2,col="gray")
  }
```

## Data for Regression Analysis

## Overview of Steps in Regression Analysis

## Estimation of Regression Function

### Example (page 15)

```{r, fig.align='center', fig.height=7,fig.width=7,message=FALSE,warning=FALSE,comment=""}

Age      <- c(20,55,30)
Attempts <- c(5, 12, 10)

fig1.9 <- data.frame(Age=Age,Attempts=Attempts)

ggplot(data= fig1.9, aes(x=Age,y=Attempts))+
	geom_point(size=4,color='gray42')+
	xlim(c(10,60))+
	ylim(c(0,15))+
  theme_bw()+
	geom_hline(yintercept = 9)+
	geom_segment(x = fig1.9[1,1], y = fig1.9[1,2], xend = fig1.9[1,1], yend = 9)+
  geom_segment(x = fig1.9[2,1], y = fig1.9[2,2], xend = fig1.9[2,1], yend = 9)+
  geom_segment(x = fig1.9[3,1], y = fig1.9[3,2], xend = fig1.9[3,1], yend = 9)+
	xlab('Age')+
	ylab('Attempts')

# Sum of Squared Deviations 

sum((Attempts - 9)^2)
```



```{r, fig.align='center', fig.height=7,fig.width=7,message=FALSE,warning=FALSE,comment=''}

mod <- lm(Attempts ~ 1 + Age,d=fig1.9)

ggplot(data= fig1.9, aes(x=Age,y=Attempts))+
	geom_point(size=4,color='gray42')+
  geom_abline(intercept = coef(mod)[1],slope=coef(mod)[2])+
	xlim(c(10,60))+
	ylim(c(0,15))+
  theme_bw()+
	geom_segment(x = fig1.9[1,1], y = fig1.9[1,2], xend = fig1.9[1,1], 
	             yend = predict(mod, newdata = data.frame(Age=c(fig1.9[1,1]))))+
	geom_segment(x = fig1.9[2,1], y = fig1.9[2,2], xend = fig1.9[2,1], 
	             yend = predict(mod, newdata = data.frame(Age=c(fig1.9[2,1]))))+
	geom_segment(x = fig1.9[3,1], y = fig1.9[3,2], xend = fig1.9[3,1], 
	             yend = predict(mod, newdata = data.frame(Age=c(fig1.9[3,1]))))+
	xlab('Age')+
	ylab('Attempts')

# Model predictions

predict(mod)

# Model Residuals

resid(mod)

# Sum of Squared Deviations 

sum((Attempts - predict(mod))^2)

sum(resid(mod)^2)

```

### Example (page 19)

```{r, fig.align='center', fig.height=5,fig.width=7,message=FALSE,warning=FALSE,comment=''}

table1.1 <- read.table(here('/data/CH01TA01.txt/'),header=FALSE)

colnames(table1.1) <- c('lot.size','work.hours')

table1.1
```

```{r, fig.align='center', fig.height=5,fig.width=7,message=FALSE,warning=FALSE,comment=''}

ggplot(table1.1,aes(x=lot.size,y=work.hours))+
  geom_point(size=3,col='gray42')+
  geom_smooth(method=lm,se=FALSE,lty=2,lwd=0.5,col='black')+
  xlim(0,150)+
  ylim(0,600)+
  xlab('Lot Size')+
  ylab('Hours')+
  theme_bw()
```


```{r, fig.align='center', fig.height=5,fig.width=7,message=FALSE,warning=FALSE,comment=''}
require(lm.beta)
require(car)

mod <- lm(work.hours ~ 1 + lot.size,d=table1.1)

Anova(mod,type=3)

summary(mod)

coef(mod)

```


```{r, fig.align='center', fig.height=5,fig.width=7,message=FALSE,warning=FALSE,comment=''}

table1.1$x_xbar <- table1.1$lot.size - mean(table1.1$lot.size)
table1.1$y_ybar <- table1.1$work.hours - mean(table1.1$work.hours)
table1.1$x_xbar_y_ybar <- table1.1$x_xbar*table1.1$y_ybar
table1.1$x_xbar.sq <- table1.1$x_xbar^2
table1.1$y_ybar.sq <- table1.1$y_ybar^2

table1.1
```

```{r, fig.align='center', fig.height=5,fig.width=7,message=FALSE,warning=FALSE,comment=''}

b1 = sum(table1.1$x_xbar_y_ybar)/sum(table1.1$x_xbar.sq)
b1

b0 = mean(table1.1$work.hours) - b1*mean(table1.1$lot.size)
b0
```


### Example (page 21)

$$ \hat{Y} = 62.37 + 35702*X$$

```{r, fig.align='center', fig.height=5,fig.width=7,message=FALSE,warning=FALSE,comment=''}


predict(mod, newdata = data.frame(lot.size=65))

```

### Table 1.2 (page 22)


```{r, fig.align='center', fig.height=5,fig.width=7,message=FALSE,warning=FALSE,comment=''}


table1.1 <- read.table(here('/data/CH01TA01.txt/'),header=FALSE)
colnames(table1.1) <- c('lot.size','work.hours')

mod <- lm(work.hours ~ 1 + lot.size,d=table1.1)

table1.1$predicted <- predict(mod)
table1.1$residuals <- resid(mod)
table1.1$residuals.squared <- resid(mod)^2

table1.1

```


### Alternative Model with Mean Centering (page 22)


```{r, fig.align='center', fig.height=5,fig.width=7,message=FALSE,warning=FALSE,comment=''}


table1.1 <- read.table(here('/data/CH01TA01.txt/'),header=FALSE)
colnames(table1.1) <- c('lot.size','work.hours')

table1.1$lot.size_centered <- table1.1$lot.size - mean(table1.1$lot.size)

mod <- lm(work.hours ~ 1 + lot.size_centered,d=table1.1)

Anova(mod,type=3)

summary(mod)

coef(mod)
```

## Estimation of Error Terms Variance


```{r, fig.align='center', fig.height=5,fig.width=7,message=FALSE,warning=FALSE,comment=''}

table1.1 <- read.table(here('/data/CH01TA01.txt/'),header=FALSE)
colnames(table1.1) <- c('lot.size','work.hours')

mod <- lm(work.hours ~ 1 + lot.size,d=table1.1)

table1.1$predicted <- predict(mod)
table1.1$residuals <- resid(mod)
table1.1$residuals.squared <- resid(mod)^2

sse = sum(table1.1$residuals.squared)
sse

mse = sum(table1.1$residuals.squared)/(nrow(table1.1)-2)
mse 

sqrt(mse)


```

## Normal Error Regression Model

### Least Square Estimation

```{r, fig.align='center', fig.height=5,fig.width=7,message=FALSE,warning=FALSE,comment='',eval=FALSE}

table1.1 <- read.table(here('/data/CH01TA01.txt/'),header=FALSE)
colnames(table1.1) <- c('lot.size','work.hours')

beta0 <- seq(50,75,.1)
beta1 <- seq(1,6,.01)

ols <- expand.grid(beta0,beta1)
colnames(ols) <- c('beta0','beta1')
ols$ssr <- NA

for(i in 1:nrow(ols)){
	p = ols[i,1] + ols[i,2]*table1.1$lot.size
	ols[i,3] = sum((table1.1$work.hours - p)^2)
}

require(lattice)

wireframe(ssr ~ beta0 * beta1,
          data = ols,
          shade=TRUE,
          screen = list(z = 40, x = -60, y=0),
	    scales = list(arrows=FALSE),
	    xlab = expression(beta[0]), 
	    ylab = expression(beta[1]), 
	    zlab = "SSR")

ols[which.min(ols$ssr),]
``` 

```{r, fig.align='center', fig.height=8,fig.width=8,message=FALSE,warning=FALSE,comment='',eval=TRUE, echo=FALSE}
ols <- read.csv(here('/data/ols.csv/'),header=TRUE)
ols[which.min(ols$ssr),]
``` 

```{r, fig.align='center', fig.height=8,fig.width=8,message=FALSE,warning=FALSE,comment='',eval=TRUE, echo=FALSE}
require(lattice)

wireframe(ssr ~ beta0 * beta1,
          data = ols,
          shade=TRUE,
          screen = list(z = 40, x = -60, y=0),
	    scales = list(arrows=FALSE),
	    xlab = expression(beta[0]), 
	    ylab = expression(beta[1]), 
	    zlab = "SSR")

``` 

### Maximum Likelihood Estimation

```{r, fig.align='center', fig.height=5,fig.width=7,message=FALSE,warning=FALSE,comment='',eval=FALSE}

table1.1 <- read.table(here('/data/CH01TA01.txt/'),header=FALSE)
colnames(table1.1) <- c('lot.size','work.hours')

beta0 <- seq(50,75,.1)
beta1 <- seq(1,6,.01)

mle <- expand.grid(beta0,beta1)
colnames(mle) <- c('beta0','beta1')
mle$loglikelihood <- NA

mse = 2383 # assumed to be known

for(i in 1:nrow(mle)){
	p = ols[i,1] + ols[i,2]*table1.1$lot.size
	mle[i,3] = 	sum(log(dnorm((table1.1$work.hours - p)/sqrt(mse))))
}

wireframe(loglikelihood ~ beta0 * beta1,
          data = mle,
          shade=TRUE,
          screen = list(z = 40, x = -60, y=0),
          scales = list(arrows=FALSE),
          xlab = expression(beta[0]), 
          ylab = expression(beta[1]), 
          zlab = "Loglikelihood")

mle[which.max(mle$loglikelihood),]

``` 

```{r, fig.align='center', fig.height=8,fig.width=8,message=FALSE,warning=FALSE,comment='',eval=TRUE, echo=FALSE}
mle <- read.csv(here('/data/mle.csv/'),header=TRUE)
mle[which.max(mle$loglikelihood),]
``` 

```{r, fig.align='center', fig.height=8,fig.width=8,message=FALSE,warning=FALSE,comment='',eval=TRUE, echo=FALSE}
require(lattice)

wireframe(loglikelihood ~ beta0 * beta1,
          data = mle,
          shade=TRUE,
          screen = list(z = 40, x = -60, y=0),
          scales = list(arrows=FALSE),
          xlab = expression(beta[0]), 
          ylab = expression(beta[1]), 
          zlab = "Loglikelihood")
``` 



## Problems

### 1.1 

No.

### 1.2

functional

### 1.3

?

### 1.4

Error component in statistical models

### 1.5

No. The expected value of Y should be the average, so there shouldn't be an error part.

$$E(Y_i) = \beta_0 + \beta_1X_i$$

### 1.6

$$Y_i = 200 + 5X_i + \epsilon_i, \epsilon_i \sim N(0,4)$$
$\beta_0 = 200$ is the expected value of Y when X = 0,
$\beta_1 = 5$ is the expected change in the value of Y for one unit increase in X


```{r, fig.align='center', fig.height=7,fig.width=7,message=FALSE,warning=FALSE}

b0 = 200
b1 = 5
        
x = 0:70
        
mean = b0 + b1*x
    
err.sd = 4
        
  plot(x,mean,type="l",
       ylim=c(150,450),
       xlim=c(0,50),
       cex=1,
       pch=19,
       xlab="X",
       ylab="Y",
       xaxt='n')
        
  axis(side=1,at=c(10,20,40))
  abline(b0,b1)
        
  dens = dnorm(seq(-3,3,.01),0,1)

  for(i in c(11,21,41)){
            
    x. = x[i] - 5*dens
    y. = mean[i]+seq(-3,3,.01)*err.sd
    points(x.,y.,type="l",lty=2)
    abline(v=x[i],lty=2,col="gray")
  }
```

### 1.7

$$Y_i = 100 + 20X_i + \epsilon_i, \epsilon_i \sim N(0,5)$$

For $X=5$, $\mu_Y = 100 + 100*5 = 200$ and $\sigma_Y = 5$. Assuming the errors are normally distributed, the probability that Y will be between 195 and 205 is about 68%. 

### 1.8

No, the expected value for the new observations would be still 104 unless you re-fit the model with the addition of this new observation.

### 1.9 

?

### 1.10

Need to see the exact coefficients for an accurate interpretation

### 1.11

No, slope only shows the relationship between the production output before and after training. The intercept difference (20) reflects the training effect. So, training increased the production output by 20.

### 1.12

a. observational, no information that the senior citizens are randomly assigned to different levels of physical activity

b. not accurate, can't make a causal statement

c. ??

d. random assignment to different levels of physical activity

### 1.13

a. observational

b. not accurate, can't make a causal statement

c. ??

d. random assignment to different levels of physical activity

### 1.14


```{r, fig.align='center', fig.height=7,fig.width=7,message=FALSE,warning=FALSE,comment=""}

set.seed(1234)

a <- data.frame(cbind(1:16,runif(16,0,1)))

a <- a[order(a[,2]),]

a$gr <- c(1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4)

a
```



```{r, fig.align='center', fig.height=7,fig.width=7,message=FALSE,warning=FALSE,eval=FALSE,echo=TRUE,comment=""}

# Simulation code to demonstrate if the above procedure works as expected

gr <- vector('list',100000)

for(i in 1:100000){
  a <- data.frame(cbind(1:16,runif(16,0,1)))
  a <- a[order(a[,2]),]
  a$gr <- c(1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4)
  gr[[i]] = a
}

assg <- c()

k <- 4

for(i in 1:100000){
  assg[i] <- gr[[i]][which(gr[[i]][,1]==k),3]
}

table(assg)

```


### 1.15

```{r, fig.align='center', fig.height=7,fig.width=7,message=FALSE,warning=FALSE,comment=""}

set.seed(1234)

a <- data.frame(cbind(1:20,runif(20,0,1)))

a <- a[order(a[,2]),]

a$gr <- c(1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4,5,5,5,5)

a
```

### 1.16

No. Least squres method doesn't care about the distribution of Y. It is still BLUE. No assumption is needed until one decides to make an inference

### 1.17

yes. why not?

### 1.18

It is not necessary that the errors will sum to zero at the population. Even the average error in population is different than zero, the residuals from a regression at the sample level can be equal to zero. See below.

```{r, fig.align='center', fig.height=7,fig.width=7,message=FALSE,warning=FALSE,comment=""}
set.seed(1234)

# Population

  b0 <- 200
  b1 <- 5
  sigma <- 5

  x   <- rnorm(10000,0,1)
  err <- rnorm(10000,50,5)
  y   <- b0+b1+err

# A random sample with size 200
  
s <- sample(1:10000,200)

s.x <- x[s]
s.y <- y[s]

fit <- lm(s.y ~ s.x)

summary(fit)

mean(resid(fit))
```

### 1.19

```{r, fig.align='center', fig.height=7,fig.width=7,message=FALSE,warning=FALSE,comment=""}
require(here)

gpa <- read.table(here('data/CH01PR19.txt'),
                  header=FALSE,
                  col.names = c('gpa','act'))

str(gpa)

head(gpa)

fit <- lm(gpa ~ act, data=gpa)

coef(fit)

summary(fit)

require(ggplot2)

ggplot(data = gpa, aes(x=act,y=gpa))+
  geom_point()+
  geom_function(fun = function(x) coef(fit)[1] + coef(fit)[2]*x)+
  theme_bw()

predict(fit)

resid(fit)

mean(resid(fit))

var(resid(fit))

sd(resid(fit))

predict(fit,
        newdata = data.frame(act = 30))

```

### 1.20

```{r, fig.align='center', fig.height=7,fig.width=7,message=FALSE,warning=FALSE,comment=""}

copier <- read.table(here('data/CH01PR20.txt'),
                  header=FALSE,
                  col.names = c('y','x'))

fit <- lm(y ~ x, data=copier)

coef(fit)

summary(fit)

ggplot(data = copier, aes(x=x,y=y))+
  geom_point()+
  geom_function(fun = function(x) coef(fit)[1] + coef(fit)[2]*x)+
  theme_bw()

predict(fit)

resid(fit)


mean(resid(fit))

var(resid(fit))

sd(resid(fit))


predict(fit,
        newdata = data.frame(x = 5))

```

### 1.21

```{r, fig.align='center', fig.height=7,fig.width=7,message=FALSE,warning=FALSE,comment=""}

freight <- read.table(here('data/CH01PR21.txt'),
                  header=FALSE,
                  col.names = c('y','x'))

fit <- lm(y ~ x, data=freight)

coef(fit)

summary(fit)

ggplot(data = freight, aes(x=x,y=y))+
  geom_point()+
  geom_function(fun = function(x) coef(fit)[1] + coef(fit)[2]*x)+
  theme_bw()

predict(fit)

resid(fit)


mean(resid(fit))

var(resid(fit))

sd(resid(fit))


predict(fit,
        newdata = data.frame(x = 1))

```


### 1.22

```{r, fig.align='center', fig.height=7,fig.width=7,message=FALSE,warning=FALSE,comment=""}

hardness <- read.table(here('data/CH01PR22.txt'),
                  header=FALSE,
                  col.names = c('y','x'))

fit <- lm(y ~ x, data=hardness)

coef(fit)

summary(fit)

ggplot(data = hardness, aes(x=x,y=y))+
  geom_point()+
  geom_function(fun = function(x) coef(fit)[1] + coef(fit)[2]*x)+
  theme_bw()

predict(fit)

resid(fit)


mean(resid(fit))

var(resid(fit))

sd(resid(fit))


predict(fit,
        newdata = data.frame(x = 40))

```

### 1.23

Refer to 1.19

### 1.24

Refer to 1.20

### 1.25

Refer to 1.21

### 1.26

Refer to 1.12

### 1.27


```{r, fig.align='center', fig.height=7,fig.width=7,message=FALSE,warning=FALSE,comment=""}

muscle <- read.table(here('data/CH01PR27.txt'),
                  header=FALSE,
                  col.names = c('y','x'))

fit <- lm(y ~ x, data=muscle)

coef(fit)

summary(fit)

ggplot(data = muscle, aes(x=x,y=y))+
  geom_point()+
  geom_function(fun = function(x) coef(fit)[1] + coef(fit)[2]*x)+
  theme_bw()

predict(fit)

resid(fit)


mean(resid(fit))

var(resid(fit))

sd(resid(fit))


predict(fit,
        newdata = data.frame(x = 40))

```


### 1.28


```{r, fig.align='center', fig.height=7,fig.width=7,message=FALSE,warning=FALSE,comment=""}

crime <- read.table(here('data/CH01PR28.txt'),
                  header=FALSE,
                  col.names = c('crime','edu'))

fit <- lm(crime ~ edu, data=crime)

coef(fit)

summary(fit)

ggplot(data = crime, aes(x=edu,y=crime))+
  geom_point()+
  geom_function(fun = function(x) coef(fit)[1] + coef(fit)[2]*x)+
  theme_bw()

predict(fit)

resid(fit)


mean(resid(fit))

var(resid(fit))

sd(resid(fit))


predict(fit,
        newdata = data.frame(edu = 80))

```

## Exercises