################################################################################
#
# Beispiel: Analyse von Tiefbohrdaten
#
# Vorlesung Lineare Modelle (Helmut Küchenhoff) 10.6.2015
# Thema: Variablenselektion
#
#       Daten zu Mineralgehalt bei Tiefbohrungen
#       file: tiefb.txt 
#
# Inhalt  : (1) Einlesen der Daten
#  	        (2) Modellbildung und Kollinearitätsanalyse
#           (3) Variablenselektion
#
# last.modified: HK 17.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")
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)

#### (3) Variablenselektion ####################################################

# volles Modell

# Anfangsmodell
reg0 <- lm(LN_CATR ~ 1,data=tiefb)

# schrittweise Variablenselection mittels AIC
stepAIC(reg0,scope=list(upper=reg, lower = reg0 ),direction = c("both"))

# schrittweise Variablenselection mittels BIC
stepAIC(reg0,scope=list(upper=reg, lower = reg0 ),direction = c("both"), 
  k = log(nrow(tiefb)))

# rückwärts Variablenselection mittels AIC
stepAIC(reg,direction = c("backward"))

# Model selection by exhaustive search, forward or backward stepwise, or sequential replacement 
library(leaps)
a <- regsubsets(LN_CATR ~ . - DEPTH - S - MGO - V  - FE2O3 - POTA- QRZ -RB,
  nbest = 1, nvmax = 10, data = tiefb)

sel <- summary(a)
coefficients(a,1:10)
sel$cp
sel$rsq
sel$bic
plot(1:10,sel$bic)
aic<-sel$bic + (2:11)*(2-log(nrow(tiefb)))
plot(1:10,aic)

### Berchnung im Modell 4 
mod4 <- lm(LN_CATR ~ C + H2O + CR + TH_COND, data = tiefb)

## Konditionszahl akzeptabel 
sqrt(kappa(mod4))
anova(mod4)
sse<-sum(mod4$residuals^2)
### aic

extractAIC(mod4,k=log(n))

### Selektion mit mallows cp 
sel$cp
plot(1:10,sel$cp)
plot(1:10,sel$cp, ylim = c(0,15))
abline(a=0,b=1)
### Der Wert bei 8 Einflussgrößen liegt erstmal unter der Hauptdiagonalen