---
title: "Jak naprawić popsutą zabawkę"
author: "Piotr J. Sobczyk"
subtitle: Transformacje zmiennych w modelach liniowych
output:
  pdf_document: default
---

```{r biblioteki, echo=FALSE, message=FALSE}
library(dplyr)
library(ggplot2)
library(lmtest)
```

> Data analysis is an artful science! It involves making subjective decisions using very objective tools!
> 
> Znalezione w notatkach wykładu Stat 501 z PennState

Na dzisiejszych zajęciach poznamy kilka sposóbów na radzenie sobie z danymi, które nie pasują do założeń
modelu liniowego. Być może między zmiennymi objaśniającymi, a zmienną objaśniającą jest zależność, ale nie jest
ona liniowa. W takim wypadku należy dokonać transformacji danych

### Transformacje zmiennych objasniających

Przede wszystkim dodajemy nowe zmienne, np. $x^2$. Mogą być też bardziej skomplikowane
przekształcenia i dyskretyzacja.

Weźmy dane dotyczące uczenia się języka obcego. Od momentu nauki sprawdzamy jaki procent słówek
pozostaje w pamięci. Dostajemy następujące dane. Interesuje nas jaka jest zależność, chcielibyśmy potrafić
przewidywać procent pamiętanych słów po danym czasie. Notabene na podobnej idei są oparte systemy takie jak
supermemo.

```{r}
wordrecall=read.csv("datasets/wordrecall.txt", sep='\t')
str(wordrecall)
```

Dopasujemy model liniowy:

```{r}
recall.lm=lm(prop~time, data=wordrecall)
summary(recall.lm)
```

Niezłe $R^2$, istotne współczynniki. A co z diagnostyką?

```{r}
par(mfrow=c(2,2))
plot(recall.lm)
par(mfrow=c(1,1))
```

Residua są normalne. Mamy natomiast problem z kiepskim dopasowaniem i niejednorodną wariancją. Mamy jedną obserwację wpływową, ale być może jest ona wynikiem złego dopasowania, więc póki co nie będziemy się nią zajmować. 

```{r}
ggplot(wordrecall, aes(x=time, y=prop))+
  geom_point() + stat_smooth(method = "lm", se = F)
```

Przeszktałcamy zmienną x, ponieważ zależność nie wygląda na liniową. Na taką potrzebę wskazuje też wykres Fitted vs Residuals. Jedna z możliwych strategii to symetryzacja zmiennej **time**. 

```{r}
MASS::boxcox(time~1, data=wordrecall)
```

```{r}
recall.lm2=lm(prop~log(time), data=wordrecall)
summary(recall.lm2)
```
Wielka różnica! Teraz $R^2$ jest bardzo wysokie! Sprawdźmy jeszcze czy naszej regresji nie zepsuliśmy.

```{r}
par(mfrow=c(2,2))
plot(recall.lm2)
par(mfrow=c(1,1))
```

Lekko niepokoić może nas obserwacja 13 z duża miarą Cooka.

Czy moglibyśmy dokonać transformacji zmiennej objaśnianej?

```{r}
MASS::boxcox(prop~time, data=wordrecall) -> boxcox.result
lambda=boxcox.result$x[which.max(boxcox.result$y)]
recall.lm3=lm(I((prop^lambda-1)/lambda)~log(time), data=wordrecall)
summary(recall.lm3)
```
Jest lepiej jeśli chodzi o $R^2$, a jak wyglądają wykresy diagnostyczne?

```{r}
qqnorm(y=MASS::stdres(recall.lm3), ylim = c(-3,3))
abline(a=0, b=1, add=TRUE)
```

Residua nie są normalne. Trochę możemy się tego spodziewać, skoro przed transformacją były normalne.

### Transformacje zmiennej objaśnianej

Uwaga. Nie robimy tego po to, żeby zmienna objaśniana miała rozkład normalny. To nie jest nasze założenie!
Założenie dotyczy rozkładu reszt.

Przekształcenie boxcoxa.

$$y' = \begin{cases} \frac{y^{\lambda}-1}{\lambda} &\mbox{dla } \lambda \neq 0 \\ 
log(\lambda) & \mbox{dla } \lambda = 0 \end{cases}$$

Sprawdźmy jak działają przekształcenia na danych dotyczących czasu trwania ciąży w ssaków.

```{r}
# codepages <- setNames(iconvlist(), iconvlist())
# x <- lapply(codepages, function(enc) try(read.csv("datasets/mammgest.txt",
#                    fileEncoding=enc,
#                    nrows=3, header=TRUE, sep="\t"), silent=T))
# unique(do.call(rbind, sapply(x, dim)))
# maybe_ok <- sapply(x, function(x) isTRUE(all.equal(dim(x), c(3,3))))
# codepages[maybe_ok]
# x[maybe_ok]
mamgest=read.csv("datasets/mammgest.txt", sep="\t", header=T, fileEncoding = "UTF-16")
str(mamgest)
```

```{r}
gestation.lm=lm(Gestation~Birthwgt, data=mamgest)
summary(gestation.lm)
```


```{r}
ggplot(mamgest, aes(x=Birthwgt, y=Gestation))+
  geom_point() + stat_smooth(method = "lm", se = F)
```

Wydaje się, że wszystko jest ok. Współczynniki duże, $R^2$ bardzo wysokie. Ale ,,diagonostyka głupcze".


```{r}
par(mfrow=c(2,2))
plot(gestation.lm)
par(mfrow=c(1,1))
```

Być może jest mały problem z normalnością residuów. Dopasowanie jest w miarę ok, ale wariancja wydaje się być niejednorodna. Zobaczmy co mówią testy.

```{r, message=FALSE}
#normalnosc reszt
shapiro.test(MASS::stdres(gestation.lm)) 
#jednorodnosc wariancji
bptest(gestation.lm)
gqtest(gestation.lm)
```

Mamy jedną obserwację wpływową, ale być może jest ona wynikiem złego dopasowania, więc póki co nie będziemy się nią zajmować. 

```{r}
MASS::boxcox(Gestation~Birthwgt, data=mamgest) -> boxcox.result
lambda=boxcox.result$x[which.max(boxcox.result$y)]
mamgest2=cbind(mamgest, GestationTransf=(mamgest$Gestation^lambda-1)/lambda)
mamgest.lm2=lm(log(Gestation)~Birthwgt, data=mamgest)
summary(mamgest.lm2)
```
Jest lepiej jeśli chodzi o $R^2$, a jak wyglądają wykresy diagnostyczne?

```{r}
par(mfrow=c(2,2))
plot(mamgest.lm2)
par(mfrow=c(1,1))
```

$R^2$ jest niższe, ale wariancja jest teraz jednorodna. Dodatkowo pozbyliśmy się obserwacji odstających. Sprawdźmy jeszcze czy nie zepsuliśmy normalności residuów

```{r, message=FALSE}
#normalnosc reszt
shapiro.test(MASS::stdres(mamgest.lm2)) 
#jednorodnosc wariancji
bptest(mamgest.lm2)
gqtest(mamgest.lm2)
```

#### Dlaczego transformacje (logarytmiczne) pomagają?

Chcemy uchwycić zależność $y$ od $x$. W modeli liniowym zakładamy, że jest ona liniowa. Co jeśli jest inaczej?

$$y=\beta_0 e^{\beta_1 x}$$
W takim wypadku nałożenie logarytmu przekształca model multikplikatywny na liniowy

$$log(y)=\beta_0 + \beta_1 x$$
Podobnie dla danych zgodnych z zależnością potęgową mamy
$$y=\beta_0 x^{\beta_1}\\ log(y)=\beta_0 + \beta_1 log(x)$$
Zaś w przypadku zależności:
$$y=\frac{x}{\beta_0 + x \beta_1}\\ \frac1y=\beta_1 + \beta_0 \frac1x$$

Reguły kciuka dla transformacji zmiennych:

1. Stężenia - przekształcenie logarytmiczne
2. Procenty - arc sin
3. Przekształcenie powinno mieć logiczne uzasadnienie. Pytanie kontrolne: czy narysowanie wykresu z danymi
w tej skali jest przekonywające?

#### Regresja wielomianowa

```{r}
data(women)
ggplot(women, aes(x=height, y=weight)) +
  geom_point() +
  stat_smooth(method = "lm", se = F, fullrange = T)
fit = lm(weight ~ height, women)
par(mfrow=c(2,3))
plot(fit, 1:6)
par(mfrow=c(1,1))
```

Bardzo złe residua, mimo, że na pierwszy rzut oka, wszystko wygląda dobrze. Problem z normalnością residuów wskazuje na problem ze skalą zmiennej objaśnianej. Dwa punkty wpływowe - najmniejszy i największy $x$.

Być może zależność nie jest liniowa. Spróbujmy dodać do modelu wyższe rzędy.

```{r}
fit2 = lm(weight ~ height + I(height^2), women)
summary(fit2)
par(mfrow=c(2,3))
plot(fit2, 1:6)
par(mfrow=c(1,1))
```
```{r}
fit3 = lm(weight ~ height + I(height^2), women, subset = -15)
summary(fit3)
par(mfrow=c(2,3))
plot(fit3, 1:6)
par(mfrow=c(1,1))
```
Zmiany są bardzo małe. Obserwacja 15 jest wpływowa, ale nie zmienia dużo.

Czy warto dodawać 3 potęgę? Niekonicznie, grozi nam zbytnie dopasowanie się do modelu (overfitting).
Zasada dotycząca umieszczania kolejnych wielomianów od $x$ jest następująca, jeśli używamy wielomianu
o stopniu $k$, to powinniśmy też używać wielomianów o wykładnikach $<k$.

### Metody radzenia sobie z obserwacjami odstającymi i wpływowymi

Dzięki diagnostyce regresji możemy zidentyfikować dwa rodzaje problematycznych obserwacji.

##### Outliers

Pierwszy to obserwacje odstajęce (ang. outliers). Są to punkty, do których regresja się źle dopasowała.
Nie muszą być one szkodliwe dla naszej regresji.

```{r}
which(rstudent(fit2)>2)
```
Dlaczego 2?

```{r}
2*(1-pnorm(2))
```


##### Influential observations

Drugi rodzaj to obserwacje wpływowe, czyli takie, których obecność wpływa na zmianę wartości współczynników
bardziej niż pozostałe punkty. Te obserwacje też nie muszą być szkodliwe dla regresji. Wpływowość mierzymy
poprzez wartość dźwigni (leverage).

```{r}
which(hatvalues(fit2)>2*ncol(fit2$model)/nrow(fit2$model))
```

Zauważmy, że wpływowość zależy wyłącznie od macierzy $X$, a nie od $y$. Czyli jest to wartość niezależna
od dopasowania modelu liniowego.

Bardzo dużym problemem są obserwacje, które są jednocześnie odstające i wpływowe. Ich niedoposowanie psuje wtedy
cały model (patrz zestaw drugi z Anscombe quartet).

Co zrobić w przypadku problemów z konkretną obserwacją? Być może jest to niepoprawna dana. Należy sprawdzić czy nie została na przykład przypadkiem niepoprawnie wprowadzona. 


Zanim usunie się daną, trzeba się zastanowić nad innym możliwymi rozwiązaniami. Może zależność nie jest liniowa i należy dokonać transformacji zmiennych objaśniających? Może brakuje nam jakiejś zmiennej objaśniającej?

Jeśli usuwamy obserwacje wpływowe o dużej wartości $x$, to powinniśmy wyraźnie zaznaczyć, że nasz model działa tylko
na wartościach $x$ z pewnego zakresu.

```{r}
infekcje=read.csv("datasets/hospital_infct_03.txt", sep="\t", fileEncoding = "UTF-16")
ggplot(infekcje, aes(x=Stay, y=InfctRsk)) + geom_point() + geom_vline(xintercept = 15, col="red")
```

Usunięcie obserwacji musi być uzasadnione! Nie jest wystarczającym powodem fakt, że model źle się dopasowuje (czyli ,,podkręcenie" wartości $R^2$). Jeśli usunięcie punktu nie wpływa na wynik analizy, to lepiej go zostawić (oczywiście umieszczając informację o obserwacji odstającej)

> First, foremost, and finally — it's okay to use your common sense and knowledge about the situation.
> 
> from Stat 501 z PennStat