---
title: "Linear Regression"
author: "Faïcel Chamroukhi"
date: "05/02/2025"
output:
  pdf_document: 
    toc: yes
  word_document: default
  html_document: default
---


## Apparts data Paris

** Some parts of this script may not be optimized as they are done in live tutorial, only **
 
```{r data, include=TRUE}

apparts<-read.csv("https://chamroukhi.com/data/apparts_Paris.csv", header=TRUE, sep=",", dec=".")

apparts

surfaces <- apparts$surfaces
prices <- apparts$prices
```

### Basic statistics

```{r, include=TRUE, echo=TRUE}
summary(surfaces)
summary(prices)
```

### Show the sample cloud

```{r data cloud, include=TRUE, echo=TRUE}

plot(surfaces, prices, xlab= "surface in m^2",ylab = "price in K euros")
```

### Fit a Linear Regression

```{r regression, include=TRUE}
# Linear Regression

x <- surfaces
y <- prices
lr_fit = lm(y ~ x)
```


### Summary of the model

```{r summary, include=TRUE}
summary(lr_fit)
```

### The fitted parameters $\hat\beta_0$ and $\hat\beta_1$ and related statistics

```{r fitted parameters, include=TRUE, echo=TRUE}
summary(lr_fit)$coefficients
```

### Plot the fitted line

```{r fitted model, include=TRUE, echo=TRUE}
y_hat = lr_fit$fitted.values

beta_0 = lr_fit$coefficients[1]
print(beta_0)
beta_1 = lr_fit$coefficients[2]
print(beta_1)

plot(x, y, xlab= "surface in m^2",ylab = "price in K euros")
lines(x, y_hat, col="red")
#lines(x, beta_0 + beta_1*x , col="blue")
#legend("topright", legend=c("data", "fitted model"),col=c("black", "red"), lty=1:2, cex=0.8)
```

### Coefficient of determination $R^2$

```{r coeff and R2,  include=TRUE, echo=TRUE}

summary(lr_fit)$r.squared

```


### Confidence Intervals for $\hat\beta_0$ and $\hat\beta_1$


```{r, include=TRUE, echo=TRUE}
alpha = 0.05
confint(lr_fit, level = 1 - alpha)
```

%```{r, include=TRUE, echo=TRUE}
%n = length(y)
%qt(1 - alpha/2, n-2) # quantile of order alpha/2 of student t with n-2 df
%```



### Statistical testing

```{r, include=TRUE, echo=TRUE}
summary(lr_fit)$coefficients
```


### Confidence interval for the regression line


```{r CI model,  include=TRUE, echo=TRUE}

seqx <- seq(min(x),max(x),length=100)

#apparts <- data.frame(surfaces, prices)
apparts <- data.frame(x, y)


Conf_Int <- predict(lr_fit, data.frame(x = seqx), interval="confidence")[,c("lwr","upr")]
summary(Conf_Int)
plot(y~x, xlab= "surface in m^2",ylab = "price in K euros",  cex=0.8) #,xlim = c(min(x),max(x)), ylim = c(min(y),max(y))
abline(lr_fit$coefficients[1],lr_fit$coefficients[2])
#matlines(seqx, cbind(Conf_Int, Pred_Int),lty=c(2,2,3,3), col=c("red","red","blue","blue"),lwd=c(2,2))
matlines(seqx, Conf_Int, lty=c(2,2), lwd=c(3,3), col=c("red","red"))

legend("bottomright",lty=c(1,3),lwd=c(1,2), c("Regression Line", paste(toString(1- alpha),"Confidence Interval")),col=c("black","red"),cex=0.8)
#lines(lr_fit$fitted.values)
#matlines(Pred_Int)
```


### Prediction interval


```{r CI/PI model,  include=TRUE, echo=TRUE}
Pred_Int <- predict(lr_fit, data.frame(x = seqx), interval="prediction")[,c("lwr","upr")]

summary(Conf_Int)

plot(y~x, xlab= "surface in m^2",ylab = "price in K euros",  cex=0.8) 

abline(lr_fit$coefficients[1],lr_fit$coefficients[2])

matlines(seqx, cbind(Conf_Int, Pred_Int),lty=c(2,2,3,3), col=c("red","red","blue","blue"),lwd=c(2,2))

legend("bottomright",lty=c(2,3),lwd=c(2,1), c("Confidence Interval","Prediction Interval"),col=c("red","blue"),cex=0.8)

#lines(lr_fit$fitted.values)
#matlines(Pred_Int)
```

### Example, for an appart of 30$m^2$, would 120K be a good deal ?

```{r CI prediction for an x*,  include=TRUE, echo=TRUE}
x0 <- 30
predict(lr_fit, data.frame(x = x0), interval="confidence")

predict(lr_fit, data.frame(x = x0), interval="prediction")

```


### Residuals
```{r,  include=TRUE, echo=TRUE}
y_hat = lr_fit$fitted.values

residuals <-  y - y_hat

#plot(x, residuals )
par(mfrow = c(1, 2))
plot(lr_fit$fitted.values, residuals, xlab = "fitted values", ylab = "residuals")
abline(0,0, col="black")
plot(lr_fit$fitted.values, y, xlab = "fitted values", ylab = "prices")
abline(0,1, col="black")
```



## Apparts data  Marseille

 
```{r, include=TRUE}

apparts<-read.csv("https://chamroukhi.com/data/apparts_Marseille.csv", header=TRUE, sep=",", dec=".")

apparts

surfaces <- apparts$surfaces
prices <- apparts$prices
```

### Basis statistics

```{r, include=TRUE, echo=TRUE}
summary(surfaces)
summary(prices)
```

### Show the sample cloud

```{r, include=TRUE, echo=TRUE}

plot(surfaces, prices, xlab= "surface in m^2",ylab = "price in K euros")
```

### Fit a Linear Regression

```{r, include=TRUE}
# Linear Regression

x <- surfaces
y <- prices
lr_fit = lm(y ~ x)
```


### Summary of the model

```{r, include=TRUE}
summary(lr_fit)
```

### The fitted parameters $\hat\beta_0$ and $\hat\beta_1$ and related statistics

```{r, include=TRUE, echo=TRUE}
summary(lr_fit)$coefficients
```

### Plot the fitted line

```{r, include=TRUE, echo=TRUE}
y_hat = lr_fit$fitted.values

beta_0 = lr_fit$coefficients[1]
print(beta_0)
beta_1 = lr_fit$coefficients[2]
print(beta_1)

plot(x, y, xlab= "surface in m^2",ylab = "price in K euros")
lines(x, y_hat, col="red")
#lines(x, beta_0 + beta_1*x , col="blue")
#legend("topright", legend=c("data", "fitted model"),col=c("black", "red"), lty=1:2, cex=0.8)
```

### Coefficient of determination $R^2$

```{r,  include=TRUE, echo=TRUE}

summary(lr_fit)$r.squared

```



### Confidence Intervals for $\hat\beta_0$ and $\hat\beta_1$


```{r, include=TRUE, echo=TRUE}
alpha = 0.05
confint(lr_fit, level = 1 - alpha)
```

### Statistical testing

```{r, include=TRUE, echo=TRUE}
summary(lr_fit)$coefficients
```


### Confidence interval for the regression line


```{r,  include=TRUE, echo=TRUE}

seqx <- seq(min(x),max(x),length=100)

#apparts <- data.frame(surfaces, prices)
apparts <- data.frame(x, y)


Conf_Int <- predict(lr_fit, data.frame(x = seqx), interval="confidence")[,c("lwr","upr")]
summary(Conf_Int)
plot(y~x, xlab= "surface in m^2",ylab = "price in K euros",  cex=0.8) #,xlim = c(min(x),max(x)), ylim = c(min(y),max(y))
abline(lr_fit$coefficients[1],lr_fit$coefficients[2])
#matlines(seqx, cbind(Conf_Int, Pred_Int),lty=c(2,2,3,3), col=c("red","red","blue","blue"),lwd=c(2,2))
matlines(seqx, Conf_Int, lty=c(2,2), lwd=c(3,3), col=c("red","red"))

legend("bottomright",lty=c(1,3),lwd=c(1,2), c("Regression Line", paste(toString(1- alpha),"Confidence Interval")),col=c("black","red"),cex=0.8)
#lines(lr_fit$fitted.values)
#matlines(Pred_Int)
```


### Prediction interval


```{r,  include=TRUE, echo=TRUE}
Pred_Int <- predict(lr_fit, data.frame(x = seqx), interval="prediction")[,c("lwr","upr")]

summary(Conf_Int)

plot(y~x, xlab= "surface in m^2",ylab = "price in K euros",  cex=0.8) 

abline(lr_fit$coefficients[1],lr_fit$coefficients[2])

matlines(seqx, cbind(Conf_Int, Pred_Int),lty=c(2,2,3,3), col=c("red","red","blue","blue"),lwd=c(2,2))

legend("bottomright",lty=c(2,3),lwd=c(2,1), c("Confidence Interval","Prediction Interval"),col=c("red","blue"),cex=0.8)

```

### Example, for an appart of 50$m^2$, would 120K be a good deal ?

```{r,  include=TRUE, echo=TRUE}
x0 <- 50
predict(lr_fit, data.frame(x = x0), interval="confidence")

predict(lr_fit, data.frame(x = x0), interval="prediction")

```



### Residuals

```{r,  include=TRUE, echo=TRUE}
y_hat = lr_fit$fitted.values

residuals <-  y - y_hat

par(mfrow = c(1, 2))
plot(lr_fit$fitted.values, residuals, xlab = "fitted values", ylab = "residuals")
abline(0,0, col="black")
plot(lr_fit$fitted.values, y, xlab = "fitted values", ylab = "prices")
abline(0,1, col="black")
```






## Apparts data  Marseille without "luxury" apparts

```{r new model,  include=TRUE, echo=TRUE}
y_new = y[y<=300]
x_new = x[y<=300]

lr_fit_new = lm(y_new ~ x_new)
y_hat_new = lr_fit_new$fitted.values

beta_0 = lr_fit_new$coefficients[1]
beta_1 = lr_fit_new$coefficients[2]

summary(lr_fit_new)$coefficients

summary(lr_fit_new)$r.squared

confint(lr_fit_new)

plot(x_new, y_new, xlab= "surface in m^2",ylab = "price in K euros")
lines(x_new, y_hat_new, col="red")
```


```{r CI prediction for an x0,  include=TRUE, echo=TRUE}
x0 <- 50
predict(lr_fit_new, data.frame(x_new = x0), interval="confidence")

predict(lr_fit_new, data.frame(x_new = x0), interval="prediction")

```



Confidence Interval

```{r CI new model,  include=TRUE, echo=TRUE}

seqx <- seq(min(x_new),max(x_new),length=1000)

#apparts <- data.frame(surfaces, prices)
apparts <- data.frame(x_new, y_new)


Pred_Int <- predict(lr_fit_new, data.frame(x_new = seqx), interval="prediction")[,c("lwr","upr")]
Conf_Int <- predict(lr_fit_new, data.frame(x_new = seqx), interval="confidence")[,c("lwr","upr")]
summary(Conf_Int)
plot(y_new~x_new, xlab= "surface in m^2",ylab = "price in K euros",  cex=0.8) #,xlim = c(min(x_new),max(x_new)), ylim = c(min(y_new),max(y_new))
abline(lr_fit_new$coefficients[1],lr_fit_new$coefficients[2])
matlines(seqx, cbind(Conf_Int, Pred_Int),lty=c(2,2,3,3), col=c("red","red","blue","blue"),lwd=c(2,2))
legend("bottomright",lty=c(2,3),lwd=c(2,1), c("Confidence Interval","Prediction Interval"),col=c("red","blue"),cex=0.8)
#lines(lr_fit_new$fitted.values)
#matlines(Pred_Int)
```

### Residuals

```{r,  include=TRUE, echo=TRUE}
y_hat = lr_fit_new$fitted.values

residuals <-  y_new - y_hat

par(mfrow = c(1, 2))
plot(lr_fit_new$fitted.values, residuals, xlab = "fitted values", ylab = "residuals")
abline(0,0, col="black")
plot(lr_fit_new$fitted.values, y_new, xlab = "fitted values", ylab = "prices")
abline(0,1, col="black")
```



### heterogeneity in the data ?
```{r clustering,  include=TRUE, echo=TRUE}
X = matrix(c(surfaces, prices),ncol = 2)
K=2
sol = kmeans(X, K)
plot(X[,1],X[,2], col = factor(sol$cluster), xlab= "surface in m^2",ylab = "price in K euros", main = "kmeans clustering")
```

