#################################################################
# Uebung zu Schaetzen und Testen I im WiSe 2015/2016
# R-Code zur Uebung 6 am 27.11.2015
# Autor: David Ruegamer, Ludwig Bothmann
#################################################################

# Aufgabe 12

# Achtung: Negative Binomialverteilung ist in R anders parametrisiert als in der Vorlesung, 
# naemlich zum Beispiel bei dnbinom mit
# size: Anzahl der Erfolge
# x: Anzahl der Misserfolge
# prob: Erfolgswahrscheinlichkeit

# Signifikanzniveau
alpha <- 0.05
# Anzahl Erfolge
k <- 20
# Erfolgswahrscheinlichkeit
pi <- 0.13


# Bestimme das (1-alpha)-Quantil der negativen Binomialverteilung 
# mit Parametern k=20 und pi0=0.13
(anzMissErfolge95 <- qnbinom(p=1-alpha,size=k,prob=pi))
# Das heisst: Wie oft wird tritt ein Misserfolg auf, bis sich mit mindestens
# 95% Wsk. 20 Erfolge ergeben
# In Summe -- weil wir nach der Angabe die "Erfolge" mitzaehlen -- ist die Schranke demnach
(c.alpha <- anzMissErfolge95 + k)

# zur Kontrolle:
pnbinom(q=c.alpha-k,size=k,prob=pi) # wieder k wegen anderer Parametrisierung abziehen
pnbinom(q=c.alpha-k,size=k,prob=pi,lower.tail=F) # s.o.
# Dass heisst: die Wahrscheinlichkeit, kleiner gleich c.alpha-k Misserfolge
# bei k=20 Erfolgen zu haben, ist pnbinom(q=c.alpha-k,size=k,prob=pi)


# Bestimmung von gamma:
(gamma <- (alpha - pnbinom(q=c.alpha-k,size=k,prob=pi,lower.tail=F))/
  dnbinom(x=c.alpha-k,size=k,prob=pi))

# Berechnung von E(Phi(X)) -> muss 0.05 ergeben
(Ephi = gamma * dnbinom(x=c.alpha-k,size=k,prob=pi) + 
  1 * pnbinom(q=c.alpha-k,size=k,prob=pi,lower.tail=F) # + 0*pnbinom(q=c.alpha-k,size=k,prob=pi) 
) # = 0.05


# Aufgabe 13
# a)

set.seed(42)
repl = 1010
trueInfl = 10
n = 50
alpha = 0.05
effectSize = 3

# Simuliere 1010x 50-dim Vektor
simMat <- matrix(rnorm(repl*n),ncol=repl) 
simMat[,1:trueInfl] <- simMat[,1:trueInfl] + effectSize

pvalsAufgabeA <- sapply((trueInfl+1):repl, function(i) t.test(simMat[,i])$p.value )
pvalsAufgabeB <- sapply(1:trueInfl,function(i) t.test(simMat[,i])$p.value )

par(mfrow=c(1,2))
hist(pvalsAufgabeA, breaks=20, xlab="") # p-Werte sollten unter H0 uniform sein
qqplot(c(1:repl)/repl,pvalsAufgabeA, xlab="theoretische Quantile", ylab="empirische Quantile")
lines(c(0,1),c(0,1), lty = "dashed", col = "blue", lwd=1.5) # unnoetig

# b)

library(parallel)
nrSims <- 100

res <- mclapply(1:nrSims, function(i){
  
  set.seed(i)
  
  # DGP
  simMat <- matrix(rnorm(repl*n),ncol=repl) 
  simMat[,1:trueInfl] <- simMat[,1:trueInfl] + effectSize
  
  pvalsAufgabeA <- sapply((trueInfl+1):repl, function(i) t.test(simMat[,i])$p.value )
  pvalsAufgabeB <- sapply(1:trueInfl,function(i) t.test(simMat[,i])$p.value )
  
  R <- sum(c(pvalsAufgabeB,pvalsAufgabeA)<=0.05)
  V <- sum(pvalsAufgabeA<=0.05)
  Q = V/R
  
  return(list(Q = Q, pvalsAufgabeA = pvalsAufgabeA, pvalsAufgabeB = pvalsAufgabeB))
  
}#, mc.cores=30  # falls mehr als ein Core vorhanden 
)

# Schaetze Erwartungswert von Q ueber Mittelwert
(FDR = mean(sapply(res, "[[", "Q")))
# Verteilung der Q
boxplot(sapply(res, "[[", "Q"))

# c)

# Adjustiere alle p-Werte nach Bonferroni

newQs <- sapply(res, function(r){
  
  newR <- sum( c(r$pvalsAufgabeB, r$pvalsAufgabeA) <=0.05 / repl)
  newV <- sum(r$pvalsAufgabeA<=0.05 / repl)
  newQ <- newV / newR
  return(newQ)
  
})

# FDR nach Korrektur
(newFDR = mean(newQs))
# Verteilung der Q
boxplot(newQs)