#########################
###### Aufgabe 1 ########
#########################

atem <- read.table("http://www.stat.uni-muenchen.de/service/datenarchiv/atem/atemwege.asc",header=T,na.strings="-1")
summary(atem)

#a)
atem.a <- atem[complete.cases(atem$pef,atem$fef50,atem$fef75,atem$sex),]
summary(atem.a)
attach(atem.a)
### Originalvariablen genügen nicht der Normalverteilung (grafisch):
plot(fef50,fef75)
# die marginalen Verteilung erscheinen eher linkssteil, daher:

### Logarithmierte Variablen:
log.fef50<-log(fef50)
log.fef75<-log(fef75)
log.fef<-cbind(log.fef50,log.fef75)
plot(log.fef50,log.fef75)

summary(log.fef)
cor(log.fef)    # Beträchtliche Korrelation (0.8846577) zwischen den Variablen

# Teil i)
#### Univariate T-tests:
t.test(log.fef50,mu=log(3),conf.level = 0.99)    # p-value = 0.08669
t.test(log.fef75,mu=log(1.4),conf.level = 0.99)  # p-value = 0.5334

#### Multivariater Test für Erwartungswert im Ein-Stichproben-Fall
#### (Kovarianz unbekannt):
n <- dim(atem.a)[1]
p <- 2
x.quer <- c(mean(log.fef50),mean(log.fef75))    # Erwartungswertvektor
mu0 <- c(log(3),log(1.4))
S <- var(log.fef)                               # geschätzte Kovarianzmatrix

# Teststatistik:
T.s <- n*(n-p)/(p*(n-1))*(x.quer-mu0)%*%solve(S) %*%(x.quer-mu0)
print(T.s)              # 11.99079
print(qf(0.99,p,n-p))   #  4.621201   => 11.99079 > 4.621201
# => Mittelwertsvektor signifikant verschieden von mu0
# p-Wert:
print(pf(T.s,p,n-p,lower.tail=F))  # 6.902e-06,lower.tail=F bedetet P[X > x]

# Man erkennt: deutlich geringerer p-Wert des multivariaten Tests als bei
# univariaten Tests (hier wird Information über starke Korrelation ausgenutzt)


# Teil ii)
#### Univariate T-tests:
T1 <- (mean(log.fef50,na.rm=T)-log(3))*sqrt(n)/sqrt(0.09)
print(abs(T1))       # 1.744487   # Betrag der Teststatistik
print(qt(0.99,n-1))  # 2.329161   # Quantil der t-Verteilung
print(pt(T1,n-1,lower.tail=F))    # p-Wert=0.9593471

T2 <- (mean(log.fef75,na.rm=T)-log(1.4))*sqrt(n)/sqrt(0.12)
print(abs(T2))       # 0.6264178  # Betrag der Teststatistik
print(qt(0.99,n-1))  # 2.329161   # Quantil der t-Verteilung
print(pt(T2,n-1,lower.tail=F))    # p-Wert=0.2655743

#### Multivariater Test für Erwartungswert im Ein-Stichproben-Fall
####(Kovarianz bekannt):

# n,p,x.quer,mu0 wie bei i)

Sigma <- matrix(c(0.09,0.07,0.07,0.12),2,2)

Z.s<-n*(x.quer-mu0)%*%solve(Sigma) %*%(x.quer-mu0)
print(Z.s)               # 8.983716
print(qchisq(0.99,p))    # 9.21034
print(pchisq(Z.s,p,lower.tail=F))  # 0.01119981
# => Mittelwertsvektor nicht verschieden von mu0


# b)

# Funktion zur Berechnung der Mahalanobis-Distanz D^2 (Ellipse)
d2.function<-function(mu1,mu2,S,xquer){

  Sinv <- solve(S)
  res <-  Sinv[1,1]*(xquer[1]-mu1)^2+
           (Sinv[1,2]+ Sinv[2,1])*(xquer[2]-mu2)*(xquer[1]-mu1)+
           Sinv[2,2]*(xquer[2]-mu2)^2
  return(res)
}

# Funktion zum Plotten der Konfidenzbereiche

draw.confi<-function(S=diag(rep(1,2)),n,alpha,xquer=rep(0,2),bonf=TRUE,...){

  p <- 2
  if (bonf){
  # Bonferroni: xi.quer+-t(n,1-alpha_i/2)*sqrt(Sii/n)
  x <- qt(1-alpha/(p*2),n-1)*sqrt(S[1,1]/n)
  y <- qt(1-alpha/(p*2),n-1)*sqrt(S[2,2]/n)
  }
  else
  {
  x <- sqrt(p*(n-1)*qf(1-alpha/(p*2),p,n-p)/(n-p))*sqrt(S[1,1]^2/n)
  y <- sqrt(p*(n-1)*qf(1-alpha/(p*2),p,n-p)/(n-p))*sqrt(S[2,2]^2/n)
  }
  # Bonferroni-Rechteck:
  plot(x=c(xquer[1]+x,xquer[1]+x,xquer[1]-x,xquer[1]-x,xquer[1]+x),
           y=c(xquer[2]+y,xquer[2]-y,xquer[2]-y,xquer[2]+y,xquer[2]+y),
           type="l",xlab="x",ylab="y",xlim=c(xquer[1]-1.5*x,xquer[1]+1.5*x),
           ylim=c(xquer[2]-1.5*y,xquer[2]+1.5*y),...)

  #Erstellen eines Grids für die Konfidenzellipsoide
  x.seq<-seq((xquer[1]-1.5*x),xquer[1]+1.5*x,length=100)
  y.seq<-seq((xquer[2]-1.5*y),(xquer[2]+1.5*y),length=100)

  #Berechnung des äußeren Produkts von x.seq und y.seq, Berechnung der
  #Mahalanobis-Distanz mit diesen Werten
  z<-outer(x.seq,y.seq,d2.function,S=S,xquer=xquer)

  #Einzeichnen der Ellipse
  ellipse<-contourLines(x.seq,y.seq,z,levels=(n-1)*p/((n-p)*n)*qf(1-alpha,p,n-p))
  lines(ellipse[[1]]$x,ellipse[[1]]$y)

}

n<-length(log.fef[complete.cases(log.fef),1])
p<-2
alpha<-0.01

draw.confi(S=S,n=n,alpha=alpha,xquer=x.quer)
points(mu0[1],mu0[2])



# c)

x3<-log(pef)
x2<-log(fef50)
x1<-log(fef75)

y<-cbind(x1-x3,x2-x3)
yquer<- c(mean(y[,1],na.rm=T),mean(y[,2],na.rm=T)) # yquer = ( -1.2580644, -0.5162401)
S <- var(y,use="complete.obs")                    
# S = 0.06934678 0.03775723
#     0.03775723 0.03254871

T.s<-(n-p)*n/(p*(n-1))*(yquer)%*%solve(S) %*%(yquer)
print(T.s)                  # 16718.75
print(qf(0.99,p-1,n-p+1))   # 6.654023
# p-Wert:
print(pf(T.s,p-1,n-p-1,lower.tail=F))     # p-Wert=0.000

# => H0 wird abgelehnt


# d)

atemvar<-cbind(fvc,pef,fef50,fef75)
summary(atemvar)
#      fvc               pef              fef50             fef75
#Min.   :  1.010   Min.   :  1.760   Min.   :  1.110   Min.   :  0.440
#1st Qu.:  1.910   1st Qu.:  4.050   1st Qu.:  2.390   1st Qu.:  1.110
#Median :  2.330   Median :  4.930   Median :  2.930   Median :  1.400
#Mean   :  2.511   Mean   :  5.161   Mean   :  3.100   Mean   :  1.498
#3rd Qu.:  3.000   3rd Qu.:  5.990   3rd Qu.:  3.643   3rd Qu.:  1.760
#Max.   :  6.050   Max.   : 12.840   Max.   :  8.760   Max.   :  5.210
#NA's   :221.000   NA's   :221.000   NA's   :221.000   NA's   :221.000

# Funktion zur Berechnung eines multivariaten Tests im Zwei-Stichproben-Fall
# für unabhängige Stichproben.

T2test <- function(varmat,gruppe,gruppenwerte=c(1,0)){

  #Argumente:
  #varmat = Matrix der Variablen
  #gruppe = Variable, nach der gruppiert werden soll
  #gruppenwerte = Ausprägung der Gruppenwerte

  #Reduktion auf vollständige Beobachtungen
  varmatgruppe<-cbind(varmat,gruppe)
  varmat<-varmat[complete.cases(varmatgruppe),]
  gruppe<-gruppe[complete.cases(varmatgruppe)]
  #Anzahl der Komponenten der Erwartungswertvektoren \mu_1 und \mu_2
  p<-dim(varmat)[2]

  #Berechnung der Mittelwertvektoren in den Gruppen
  x1.quer<-apply(varmat[gruppe==gruppenwerte[1],],MARGIN=2,FUN=mean)
  x2.quer<-apply(varmat[gruppe==gruppenwerte[2],],MARGIN=2,FUN=mean)

  #Berechnung Stichprobenumfänge in den Gruppen
  n1<-sum(gruppe==gruppenwerte[1])
  n2<-sum(gruppe==gruppenwerte[2])

  #Berechnung der Kovarianzmatrizen in den Gruppen
  S1<-var(varmat[gruppe==gruppenwerte[1],])
  S2<-var(varmat[gruppe==gruppenwerte[2],])

  #Berechnung der gepoolten Kovarianzmatrix
  S<-1/(n1+n2-2)*((n1-1)*S1+(n2-1)*S2)

  #Teststatistik T^2 (vgl. Aufgabe 1)
  T2<-n1*n2/(n1+n2)*(x1.quer-x2.quer)%*%solve(S)%*%(x1.quer-x2.quer)
  # p-Wert
  pwert<-pf((n1+n2-p-1)/((n1+n2-2)*p)*T2,p,n1+n2-p-1,lower.tail=F)
  erg<-list("T2"=T2,"pwert"=pwert)
  return(erg)
}


# Test für Gruppenunterschiede bzgl. Geschlecht (1: männlich, 2: weiblich)

summary.factor(sex)         #804 745

# Univariat:
t.test(fvc[sex==1],fvc[sex==2],var.equal=T)         # p-value = 9.978e-16
t.test(pef[sex==1],pef[sex==2],var.equal=T)         # p-value = 4.74e-09
t.test(fef50[sex==1],fef50[sex==2],var.equal=T)     # p-value = 0.01554
t.test(fef75[sex==1],fef75[sex==2],var.equal=T)     # p-value = 0.4475
# => bis auf fef75 signifikante Unterschiede zum 0.05-Niveau.

# Multivariat:
T2test(atemvar[,3:4],gruppe=sex,gruppenwerte=c(1,2))    # p-value = 0.0009796528
# => Unterschiede "hochsignifikant"


detach(atem)