############################################################
########## Úvod do korelační a regresní analýzy ############
######### Martina Litschmannová, Adéla Vrtková #############
############################################################

##### Vícenásobná lineární regrese = korelační analýza + regresní analýza
?stackloss
head(stackloss)

############# Korelační analýza ############################
# 1. Explorační analýza korelačního pole
  # korelační koeficienty
  # identifikace odlehlých pozorování
  # parciální korelační koeficienty
 
############# Regresní analýza ############################
# 1. Explorační analýza korelačního pole
  # korelační koeficienty
  # identifikace odlehlých pozorování
  # parciální korelační koeficienty
  # odhad typu regresní funkce (není-li známý)
  # multikolinearita
# 2. Odhad koeficientů regresní funkce
# 3. Verifikace modelu
  # F-test, dílčí t-testy
  # identifikace vlivných bodů
# 4. Ověření kvality modelu
# 5. Analýza reziduí
# 4. Predikce
####################################################################

#########################################################
# 1. Explorační analýza korelačního pole
#########################################################

#### korelační koeficienty, korelační matice ####

# náhled na závislost mezi dvojicemi sledovaných proměnných - tzv. Scatterplot matrix
par(mfrow = c(1,1))
pairs(stackloss) # matice korelačních polí (bodových grafů)

# funkce, která do jednoho grafického okna vykreslí krabicové grafy, histogramy, QQ-grafy
# navíc provede test normality
my_exp_reg = function(df){
  par(mfrow = c(3,4))
  boxplot(df$stack.loss,xlab = "Stack loss",horizontal = T)
  boxplot(df$Air.Flow,xlab = "Air Flow",horizontal = T)
  boxplot(df$Water.Temp,xlab = "Water Temperature",horizontal = T)
  boxplot(df$Acid.Conc.,xlab = "Acid Concentration",horizontal = T)
  hist(df$stack.loss,xlab = "Stack loss",main = "")
  hist(df$Air.Flow,xlab = "Air Flow",main = "")
  hist(df$Water.Temp,xlab = "Water Temperature",main = "")
  hist(df$Acid.Conc.,xlab = "Acid Concentration",main = "")
  qqnorm(df$stack.loss, main = "qqplot - Stack Loss")
  qqline(df$stack.loss)
  qqnorm(df$Air.Flow, main = "qqplot - Air Flow")
  qqline(df$Air.Flow)
  qqnorm(df$Water.Temp, main = "qqplot - Water Temperature")
  qqline(df$Water.Temp)
  qqnorm(df$Acid.Conc., main = "qqplot - Acid Concentration")
  qqline(df$Acid.Conc.)

r1 = paste("test normality (Stack Loss):",round(shapiro.test(df$stack.loss)$p.value,3))
r2 = paste("test normality (Air Flow):",round(shapiro.test(df$Air.Flow)$p.value,3))
r3 = paste("test normality (Water Temperature):",round(shapiro.test(df$Water.Temp)$p.value,3))
r4 = paste("test normality (Acid Concentration):",round(shapiro.test(df$Acid.Conc.)$p.value,3))

print(r1)
print(r2)
print(r3)
print(r4)
}

my_exp_reg(stackloss)

# náhled na korelační matice s Pearsonovými a Spearmanovými korel. koeficienty
cor(stackloss, method = "pearson")
cor(stackloss, method = "spearman")

# ke korelačním koeficientům potřebujeme znát i jejich významnost (p-hodnoty)
# install.packages("Hmisc")
library(Hmisc)
rcorr(as.matrix(stackloss), type = "spearman") #funkce vrátí jak korelační matici, tak p-hodnoty
korel_mat = round(rcorr(as.matrix(stackloss), type = "spearman")[[1]],3)
korel_mat_p = round(rcorr(as.matrix(stackloss), type = "spearman")[[3]],3)

korel_mat
korel_mat_p

# korelační matici společně s významností korel. koeficientů lze i graficky zobrazit
par(mfrow = c(1,1))
# install.packages("corrplot")
library(corrplot)
corrplot(korel_mat) # korelogram lze vylepšit - úprava barev, pouze horní (dolní) část, apod.

corrplot(korel_mat,
         method = "color", # způsob vizualizace (zkuste si např. "square", "ellipse")
         addCoef.col = "black",  # přidá hodnotu korel. koeficientu
         type = "upper",  # horní trojúhelník
         order = "hclust",
         tl.col="black", tl.srt=45, # Nastavení barvy a rotace textových popisků
         # kombinace grafu s významnosti korelačního koeficientu
         p.mat = korel_mat_p, sig.level = 0.05, insig = "blank", 
         diag=FALSE # skrytí korelačních koeficientů na diagonále
         )

# změna barev

corrplot(korel_mat,
         method = "circle",
         addCoef.col = "black",
         type = "upper",
         order = "hclust",
         tl.col="black", tl.srt=30, # Nastavení barvy a rotace textových popisků
         # kombinace grafu s významnosti korelačního koeficientu
         p.mat = korel_mat_p, sig.level = 0.05, insig = "blank", 
         diag=FALSE, # skrytí korelačních koeficientů na diagonále
         col = topo.colors(n = 8) # změna barev
         )

# Můžeme si napsat vlastní funkci, kde vše pak získáme najednou
# Funkce pro výpočet a vizualizaci korelační matice
# Funkce pracuje s balíčky Hmisc a corrplot, které musí být nainstalovány a načteny
kor_mat_funkce = function(df,typ = "pearson",alpha = 0.05){
  # df ... datová matice (data frame),
  # type ... "pearson" nebo "spearman"
  # alpha ... kritická hladina významnosti
  korel_mat = round(rcorr(as.matrix(df), type = typ)[[1]],3)
  korel_mat_p = round(rcorr(as.matrix(df), type = typ)[[3]],3)
  par(mfrow = c(1,1))
  corrplot(korel_mat,
           method = "circle",
           addCoef.col = "black",
           type = "upper",
           order = "hclust",
           tl.col="black", tl.srt=30, # Nastavení barvy a rotace textových popisků
           # kombinace grafu s významnosti korelačního koeficientu
           p.mat = korel_mat_p, sig.level = alpha, insig = "blank", 
           diag=FALSE, # skrytí korelačních koeficientů na diagonále
           col = topo.colors(n = 8) # změna barev
  )
  text1 = paste0("korelační matice (",typ,"):")
  text2 = paste0("korelační matice (",typ,") - p - hodnoty:")
  return(list(text1,korel_mat,text2,korel_mat_p))
}

kor_mat_funkce(stackloss,"spearman",0.05)

### Identifikace odlehlých pozorování ####

# Jak najít odlehlá pozorování v jednotlivých proměnných?
stackloss_out = stackloss # POZOR! Je vhodné pracovat s kopii původních dat!
# např. pro 1. sloupec - proměnná Air.Flow
boxplot.stats(stackloss_out[,1])
out = boxplot.stats(stackloss_out[,1])$out

# Jak odstranit odlehlá pozorování z proměnné Air.Flow? (nahradit je NA)
stackloss_out[stackloss_out[,1] %in% out,1] = NA

# Identifikace odlehlých pozorování v jednotlivých proměnných
# a jejich nahrazení NA

stackloss_out = stackloss
for (i in 1:4){
  out = boxplot.stats(stackloss[,i])$out
  stackloss_out[stackloss_out[,i] %in% out,i] = NA
  }

# Výběr kompletních záznamů - dle rozhodnutí datového analytika
stackloss_complete = na.omit(stackloss_out)

# Explorační a korelační analýza kompletních dat po odstranění NA
my_exp_reg(stackloss_complete)
kor_mat_funkce(stackloss_complete,"spearman",0.05)

# Místo vyloučení celého pozorování, lze nahradit (imputovat) chybějící hodnoty. Metod pro nahrazení je více,
# např. průměrem, mediánem, lineární aproximací, ...

# Pro zájemnce
# balíček mice (Multivariate Imputation by Chained Equations) 
# je balíček R, který poskytuje pokročilé funkce pro nahrazení chybějících
# hodnot. Využívá nepříliš známý způsob nahrazení NA ve dvou krocích 
# pomocí funkce mice() k sestavení modelu a funkce complete() k vygenerování 
# datového souboru (dále df) bez NA. Funkce mice(df) vytváří několik 
# kopií df, z nichž každá má aplikován jiný způsob nahrazení NA. 
# Funkce complete() vrátí jeden nebo několik df, přičemž výchozí je první.

##### Parciální korelační koeficienty ##########
# (Srovnejte Spearmanův korelační koeficient pro stack.loss a Air.Flow 
# s příslušným parciálním korelačním koeficientem s vyloučením vlivu Water.Temp a Acid.Conc.)

cor.test(stackloss_complete$stack.loss,stackloss_complete$Air.Flow,method = "spearman")
# install.packages("ppcor")
library(ppcor)
pcor.test(stackloss_complete$stack.loss,stackloss_complete$Air.Flow,
          c(stackloss_complete$Water.Temp,stackloss_complete$Acid.Conc.),
          method = "spearman")

###################################################################################
### Další část explorační analýzy děláme jen v případě, že chceme pokračovat regresní analýzou
###################################################################################

# Odhad typu regresní funkce
par(mfrow = c(1,3))
plot(stack.loss ~ Air.Flow,data = stackloss_complete)
plot(stack.loss ~ Water.Temp,data = stackloss_complete)
plot(stack.loss ~ Acid.Conc.,data = stackloss_complete)

# Budeme hledat model ve tvaru: 
# stack.loss = konst. + b0*Air.Flow + b1*Water.Temp + b2*Acid.Conc. + e

# Ověření multikolinearity (existence lineární závislosti mezi prediktory)
# abs. hodnota korel. koeficientu větší než 0,8 je považována za problematickou
pairs(stackloss_complete[,-4])
cor(stackloss_complete[,-4],method = "pearson") 
kor_mat_funkce(stackloss_complete[,-4],"pearson",0.05)
# Nepředpokládáme, že mezi prediktory existuje multikolinearita. 

###################################################################################
# 2. Odhad koeficientů regresní funkce - aplikace metody nejmenších čtvercůOdhad typu regresního modelu 
###################################################################################

model = lm(stack.loss ~ Air.Flow + Water.Temp + Acid.Conc., data = stackloss_complete) 
# nebo
model = lm(stack.loss ~ . , data = stackloss_complete)
model

# ve funkci lm() je parametr na.action nastaven na na.omit
# srovnejte model_out a model
model_out = lm(stack.loss ~ . , data = stackloss_out)
model_out

# Jak vypustit konstantu z modelu?
model_2 = lm(stack.loss ~ 0 + Air.Flow + Water.Temp + Acid.Conc., data = stackloss_complete)
model_2
summary(model_2)
###################################################################################
# 3. Verifikace modelu
###################################################################################

####### F-test, dílčí t-testy ######################
summary(model)

# Zjednodušení modelu dle výsledků dílčích t-testů -> vypuštění Acid.Conc
model_3 = lm(stack.loss ~ Air.Flow + Water.Temp, data = stackloss_complete)
summary(model_3)

# Model je dobře specifikován (p-hodnota<0.001, F-test).
# "Všechny použité prediktory jsou statisticky významné (viz dílčí t-testy)."

#########################################################################
# Máme:
# model - model se všemi původními regresory, s abs. členem
# model_2 - model se všemi původními regresory, bez abs. členu
# model_3 - model bez Acid.Conc proměnné, s abs. členem
##########################################################################


###### Identifikace vlivných bodů ##########
im = influence.measures(model_3)
im
summary(im)  # potenciální vlivné body

# vyloučíme pozorování "4" "21"
stackloss_complete_3_2 = stackloss_complete[rownames(stackloss_complete)!= "4" & rownames(stackloss_complete)!= "21",]
model_3_2 = lm(stack.loss ~ Air.Flow + Water.Temp , data = stackloss_complete_3_2)
summary(model_3_2)

#########################################################################
# Máme:
# model - model se všemi původními regresory, s abs. členem
# model_2 - model se všemi původními regresory, bez abs. členu
# model_3 - model bez Acid.Conc proměnné, s abs. členem
# model_3_2 - model bez Acid.Conc proměnné, s abs. členem - z dat, z nichž byly odstraněny vlivné body
##########################################################################

###################################################################################
# 4. Posouzení kvality modelu
###################################################################################

# Vizuální posouzení kvality modelu
plot(stackloss_complete_3_2$stack.loss,model_3_2$fitted.values,
     col = "red",
     pch = 16,
     xlab = "pozorované hodnoty vysvětlované proměné",
     ylab = "predikované hodnoty")
points(stackloss_complete$stack.loss,model$fitted.values, pch = 16)
lines(0:20,0:20)

# R - kvadrát - index determinace
summary(model)
summary(model_3_2)

# Akaikeho kritérium
AIC(model) # čím menší, tím lepší
AIC(model_3_2)

# Bayesovo informační kriterium 
# (vyšší penalizace pro počet nezávisle proměnných)
BIC(model) # čím menší, tím lepší
BIC(model_3_2)

# Další kritéria pro posouzení kvality modelu (MAE, RMSE, ...)
install.packages("DMwR")
library(DMwR)
regr.eval(stackloss_complete$stack.loss,model$fitted.values)
regr.eval(stackloss_complete_3_2$stack.loss,model_3_2$fitted.values)


###################################################################################
# 5. Analýza reziduí
###################################################################################

par(mfrow = c(2,2))
plot(model_3)  # POZOR! Je-li vstupem funkce plot regresní model, 
               # výstupem jsou 4 grafy! Pro nás jsou zajímavé první dva.

# ověření normality reziduí
shapiro.test(model_3$residuals) # očekávali jsme - viz grafy regr. diagnostiky

# ověření nulové střední hodnoty reziduí
t.test(model_3$residuals)

# ověření, zda rezidua nejsou autokorelovaná
install.packages("lmtest")
library(lmtest)
dwtest(model_3)

# Breuschův - Paganův test homoscedasticity reziduí
install.packages("car")
library(car)
ncvTest(model_3) 

###################################################################################
# 6. Predikce
###################################################################################
model_3_2
newdata = data.frame(Air.Flow = 60,Water.Temp = 20) # zadáme nové hodnoty vysvětlujících proměnných do datového rámce
predict(model_3_2,newdata,interval = "confidence")
predict(model_3_2,newdata,interval = "predict")