################################################################################
#
# Beispiel: Analyse von Tiefbohrdaten
#
# Vorlesung Lineare Modelle (Helmut Küchenhoff) 3.6. 2015
# Thema: Diagnose 
#
#       Daten zu Mineralgehalt bei Tiefbohrungen
#       file: tiefb.txt 
#
# Inhalt  : (1) Einlesen der Daten
#  	        (2) Kollinearitätsanalyse
#       
#
# last.modified: HK 3.6.2015
#
################################################################################

### (1) Einlesen der Daten #####################################################

# Einlesen der Daten und kurze Kontrolle
tiefb = read.table('tiefb.txt', header = TRUE)
# Variablen:
# DEPTH - Tiefe der Probe
# LN_CATR - logarithmiert
# H20, C, QRZ, ... - verschiedene weitere Mineralien in der Probe
# DENSITY - Dichte der Probe

# kurze Kontrolle
head(tiefb)
plot(LN_CATR ~ DEPTH, data=tiefb, type='l')
n <- nrow(tiefb)
#### (2) Modellbildung #########################################################

library(MASS)
library(car)

# einfache multiple Regression
reg <- lm(LN_CATR ~ H2O + C + QRZ + AL2O3 + FE2O3 + MGO + CAO + NA2O + S + SR + 
  RB + ZR + CR + V + CU + POTA + URANIUM + THORIUM + SUS + DENSITY + TH_COND, 
  data=tiefb)
# oder viel kürzer
reg <- lm(LN_CATR ~ . - DEPTH  , data = tiefb)  

# Kollinearitätsanalyse
vif(reg)
kappa(reg)
### Es fallen einige Vraibalen auf:AL2O3     FE2O3  
cor(tiefb$AL2O3,tiefb$FE2O3)

library(corrplot)
### Correlation der Einflussgrößen 
xgr<-(tiefb[,-c(1,2)])
corrplot(cor(xgr),method = "color", order = "hclust")
### Alle einflussgrößen mit cor>0.8
xgr1<-xgr[,vif(reg) > 10]
corrplot(cor(xgr1),method = "number",order = "hclust")
# entferne Kollineare Terme aus dem Modell und belasse Stellvertreter
reg = lm(LN_CATR ~ . - DEPTH - S - MGO -V - FE2O3 , data = tiefb)        
vif(reg)
### Konditionszahl 
sqrt(kappa(reg))
corrplot(cor(xgr), method="color", order="hclust")
dev.off()
reg = lm(LN_CATR ~ . - DEPTH - S - MGO -V - FE2O3 - QRZ , data=tiefb)   
# entferne weitere Kollineare Terme aus dem Modell und belasse Stellvertreter
vif(reg)
reg = lm(LN_CATR ~ . - DEPTH - S - MGO - V  - FE2O3 - POTA- QRZ -RB , data = tiefb)
vif(reg)
plot(reg)
plot(rstandard(reg)[1:386],rstandard(reg)[2:387])
plot(cook.distance(reg))
### Autokorrelation wegen Tiefe 
cor(rstandard(reg)[1:386],rstandard(reg)[2:387])
## Durbin Watson Test 
dwt(reg)