#  Copyright (C) 2011-2012 Sandra Ligges
#
#  This program is free software; you can redistribute it and/or modify it under the terms
#  of the GNU General Public License as published by the Free Software Foundation;
#  either version 2 of the License, or (at your option) any later version.
#
#  This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
#  without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
#  See the GNU General Public License for more details.
#
#  A copy of the GNU General Public License is available in the same folder as this file.


#######################################################################################################################
### Funktion zur Erstellung der Verlaufsdiagramme für die Ergebnisse der Intervallschätzverfahren
#######################################################################################################################
SimEstHR.KIs <- function(file.erg1, file.erg2) {
    
    ### Vorbereitungen
    ## Definition der betrachteten (logarithmierten) Hazard-Ratios und eines Indizes für die Szenarien
    haz.rat <- c(1/3, 1/2, 2/3, 1/1.2, 1, 1.2, 1.5, 2, 3) 
    log.haz.rat <- log(haz.rat)
    haz.rat.l <- length(haz.rat)
    idx.sz <- 1:50
    
    ## Laden des Versuchsplans
    load(file = paste(file.erg1, "/SimEstHR_Info_Vplan.RData", sep = ""))

    ## Spezifikationen, die für den Titel der Abbildungen benötig werden
    Vplan.Zens <- ifelse(Vplan[, 6] == "keine", "nein", "ja")
    Vplan.Bdg <- ifelse(Vplan[, 7] == "keine", "nein", "ja")
    Vplan.Stoer <- ifelse(Vplan[, 8] == "keiner", "nein", "ja")
    
    ## Definition von Vektoren, die für ein gegebenes Szenario und Hazard-Ratio 
    ## die Überdeckungsraten und durchschnittlichen Breiten beinhalten sollen
    cox.ud <- numeric(haz.rat.l)
    was.ud <- numeric(haz.rat.l)
    sl.ud <- numeric(haz.rat.l)
    cox.b <- numeric(haz.rat.l)
    was.b <- numeric(haz.rat.l)
    sl.b <- numeric(haz.rat.l)
    
    ## Berechnung der Überdeckungserfolgsgrenze
    runs <- 10000
    ud.grenze <-  qnorm(0.05) * sqrt(0.95 * (1 - 0.95) / runs) + 0.95 
    
    ## Folgende Szenarien der alten Nummerierung bilden jeweils das erste Szenario in den zu zeichnenden Abbildungen
    sz.reihf <- c(3, 4, 1, 2, 7, 8, 5, 6, 11, 12, 9, 10, 15, 16, 13, 14, 49) 
    
    
    ### Zeichnung der Grafiken für je zwei bis drei der Szenarien in einer getrennten Datei ###########################
    ## Setzen eines Zählers, der die Nummern der jeweils ersten Szenarios in den Abbildungen durchläuft
    z <- -2
    sz.reihf <- c(3, 4, 1, 2, 7, 8, 5, 6, 11, 12, 9, 10, 15, 16, 13, 14, 49)
    
    ## Durchlauf der einzelnen Abbildungen
    for(i in sz.reihf) {
    
        ## Erstellung der Szenarioauswahl für eine Abbildung sowie
        ## Erzeugen der Überschriften der zwei bis drei Grafiken in einer Abbildung
        if(i < 49) {
            sz.ausw <- i + c(0, 16, 32)
            main.text <- c("Exp(variierend)/Exp(0.0693) (Gruppe 0/Gruppe 1)", 
                "Weib(0.5, variierend)/Weib(0.5, 0.2192) (Gruppe 0/Gruppe 1)", 
                "Gomp(0.3, variierend)/Gomp(0.3, 0.0109) (Gruppe 0/Gruppe 1)")
        }else{ 
            sz.ausw <- c(49, 50)
            main.text <- c("Weib(0.5, variierend)/Weib(0.5, 0.2192) (Gruppe 0/Gruppe 1), Störeffekt: 2", 
                "Weib(0.5, variierend)/Weib(0.5, 0.2192) (Gruppe 0/Gruppe 1), Störeffekt: 4")
        }
        
        ## Setzen eines weiteren Zählers für die betrachteten Szenarien innerhalb einer Abbildung 
        zaehler <- 0    
        
        ## Definition einer Liste, die alle zu plottenden Werte einer Abbildung tragen soll
        L <- vector(mode = "list", length = length(sz.ausw))
        
        ## Durchlauf der zwei bzw. drei Szenarien einer Szenarioauswahl für eine Abbildung
        for(j in sz.ausw) {
            
            # Laden der Datei des entsprechenden Szenarios
            load(paste(file.erg1, "/SimEstHR_Sz", idx.sz[j], ".RData", sep = ""))
            
            # Extrahierung der Konfidenzintervallgrenzen der drei Intervalltypten aller Hazard-Ratios eines Szenarios
            # und Berechnung der zugehörigen Überdeckungsraten und durchschnittlichen Breiten der Intervalle
            zaehler <- zaehler + 1 
            for(l in seq(along = haz.rat)) {
                konf <- sim.erg[[l]][, 14:19]
        
                cox.ud[l] <- mean((konf$Cox.KI.u < log.haz.rat[l]) & (konf$Cox.KI.o > log.haz.rat[l]), na.rm = TRUE)
                was.ud[l] <- mean((konf$Wassmer.KI.u < log.haz.rat[l]) & (konf$Wassmer.KI.o > log.haz.rat[l]), 
                    na.rm = TRUE)
                sl.ud[l] <- mean((konf$SL.KI.u < log.haz.rat[l]) & (konf$SL.KI.o > log.haz.rat[l]), na.rm = TRUE)
        
                cox.b[l] <- mean(konf$Cox.KI.o - konf$Cox.KI.u, na.rm = TRUE)
                was.b[l] <- mean(konf$Wassmer.KI.o - konf$Wassmer.KI.u, na.rm = TRUE)
                sl.b[l] <- mean(konf$SL.KI.o - konf$SL.KI.u, na.rm = TRUE)
            }
            L[[zaehler]] <- list(cox.ud = cox.ud, was.ud = was.ud, sl.ud = sl.ud, 
                cox.b = cox.b, was.b = was.b, sl.b = sl.b)
        }
        
        z <- z + 3
        ## Öffnen eines pdfs
        pdf(paste(file.erg2, "/SimEstHR_KIs_", z, ".pdf", sep = ""), width = 21/2.54, height = 29.7/2.54)
        
        ## allgemeine Einstellungen für das Layout der Abbildungen
        layout(matrix(c(3, 3, 1:2, 6, 6, 4:5, 9, 9, 7:8), ncol = 2, byrow = TRUE), heights = c(0.1, 1, 0.1, 1, 0.1, 1))
        par(lwd = 1.2, cex = 0.6, cex.axis = 1.4, cex.lab = 1.4, cex.main = 1.45, 
            mar = c(4.5, 5, 0.5, 1), oma = c(0.1, 0.1, 3.5, 0.1))  
        
        ## Zeichnen der Abbildungen mit den 4-6 Grafiken
        zaehler <- 0        
        for(j in sz.ausw) {
            zaehler <- zaehler + 1  
            
            # Extrahierung der bereits berechneten zu plottenden Werte eines Szenarios aus der Szenarioauswahl
            cox.ud <- L[[zaehler]]$cox.ud
            was.ud <- L[[zaehler]]$was.ud
            sl.ud <- L[[zaehler]]$sl.ud
            cox.b <- L[[zaehler]]$cox.b
            was.b <- L[[zaehler]]$was.b
            sl.b <- L[[zaehler]]$sl.b
            
            # Berechnung des zu plottenden Bereichs der y-Achse in der aktuellen Grafik der Überdeckungsraten
            min.ud <- max(0.85, min(c(0.915, sapply(L, "[[", "cox.ud"), sapply(L, "[[", "was.ud"), 
                sapply(L, "[[", "sl.ud")), na.rm = TRUE))
            max.ud <- 0.985
            
            # Zeichnen der Grafik der Überdeckungsraten für das aktuelle Szenario
            plot(log.haz.rat, cox.ud, type = "b", ylim = c(min.ud, max.ud), xlim = c(-1.25, 1.25), xaxt = "n",
                xlab = expression(Wahres ~~ ln(omega)), ylab = "Überdeckungsrate", main = "", lwd = 1.5)
            axis(1, at = log.haz.rat, labels = round(log.haz.rat, 2))
            
            lines(log.haz.rat, was.ud, type = "b", lty = "dashed", lwd = 1.5)
            lines(log.haz.rat, sl.ud, type = "b", lty = "dotted", lwd = 1.8)
            
            # Einzeichnen einer Hilfslinie für das Konfidenzniveau und einer für die Überdeckungserfolgsgrenze
            abline(h = 0.95)
            abline(h = ud.grenze, lty = "dotted")   
            
            # Erstellen einer Legende für die Niveaugrenze und die Überdeckungserfolgsgrenze
            if((j <= 16) || j == 49) {
                legend("topleft", horiz = TRUE, legend = c("Überdeckungserfolgsgrenze .9464",
                    expression("Niveau" ~ alpha == ".95")), lty = c("dotted", "solid"), lwd = 1.2, cex = 1.1)
            } 
            
            # Punkte, die außerhalb des zu plottenden y-Bereiches liegen,
            # werden als entsprechender Text in die Grafiken eingetragen.
            min.g <- par("usr")[3]
            if(!(j %in% seq(2, 48, by = 2))) {
                if(any(cox.ud < min.g))
                    text(log.haz.rat[cox.ud < min.g], 0.900, 
                        as.expression(lapply(round(cox.ud[cox.ud < min.g], 2), function(x)
                            substitute(KI[C] * ": " * s, list(s = x)))), 
                        cex = 1.1)
                if(any(was.ud < min.g))
                    text(log.haz.rat[was.ud < min.g], 0.875, 
                        as.expression(lapply(round(was.ud[was.ud < min.g], 2), function(x)
                            substitute(KI[W] * ": " * s, list(s = x)))), 
                        cex = 1.1)
                if(any(sl.ud < min.g))
                    text(log.haz.rat[sl.ud < min.g], 0.850,
                        as.expression(lapply(round(sl.ud[sl.ud < min.g], 2), function(x)
                            substitute(KI[L] * ": " * s, list(s = x)))), 
                        cex = 1.1)
            }
            else{
                if(any(cox.ud < min.g))
                    text(log.haz.rat[cox.ud < min.g], 0.900, 
                        as.expression(lapply(round(cox.ud[cox.ud < min.g], 2), function(x)
                            substitute(KI[C]^B * ": " * s, list(s = x)))), 
                        cex = 1.1)
                if(any(was.ud < min.g))
                    text(log.haz.rat[was.ud < min.g], 0.875, 
                        as.expression(lapply(round(was.ud[was.ud < min.g], 2), function(x)
                            substitute(KI[W]^B * ": " * s, list(s = x)))), 
                        cex = 1.1)
                if(any(sl.ud < min.g))
                    text(log.haz.rat[sl.ud < min.g], 0.850,
                        as.expression(lapply(round(sl.ud[sl.ud < min.g], 2), function(x)
                            substitute(KI[L]^B * ": " * s, list(s = x)))), 
                        cex = 1.1)
            }
            
            # Berechnung des zu plottenden Bereichs der y-Achse in der aktuellen Grafik der Konfidenzintervallbreiten
            min.b <- min(c(sapply(L, "[[", "cox.b"), sapply(L, "[[", "was.b"), sapply(L, "[[", "sl.b")), na.rm = TRUE)
            max.b <- max(c(sapply(L, "[[", "cox.b"), sapply(L, "[[", "was.b"), sapply(L, "[[", "sl.b")), na.rm = TRUE)
            min.b <- min.b - 0.1 * (max.b - min.b)
            
            # Zeichnen der Grafik der durchschnittlichen Konfidenzintervallbreiten für das aktuelle Szenario
            plot(log.haz.rat, cox.b, type = "b", ylim = c(min.b, max.b), xaxt = "n",
                xlab = expression(Wahres ~~ ln(omega)), ylab = "Durchschnittliche Breite", main = "", lwd = 1.5)
            axis(1, at = log.haz.rat, labels = round(log.haz.rat, 2))
            
            lines(log.haz.rat, was.b, type = "b", lty = "dashed", lwd = 1.5)
            lines(log.haz.rat, sl.b, type = "b", lty = "dotted", lwd = 1.8)
            
            # Erstellen einer Legende für die drei abbgebildeten Intervalltypen
            if((j > 32 && j < 49) || j == 50) {
                legend("bottom", xjust = 0.5, yjust = 0, horiz = TRUE, legend = c("Cox", "Wassmer", "SL"), 
                    lty = c("solid", "dashed", "dotted"), lwd = 1.8, cex = 1.2)
            }
            
            # Beschriftung der beiden Grafiken für das aktuelle Szenario    
            op <- par(mar = c(0, 0, 2.5, 0), xpd = NA)
            plot.new()
            title(main = main.text[zaehler])
            par(op)
        }
         
        ## Erstellen eines Titels für die aktuelle Abbildung  
        if(i < 49) {
            title(paste("Stichprobenumfang: ", Vplan[i, 4], "/", Vplan[i, 5], " (Gruppe 0/Gruppe 1),", 
                "\n Zensierungen: ", Vplan.Zens[i], ", Bindungen: ", Vplan.Bdg[i], ", Störeffekt: ", 
                Vplan.Stoer[i], sep = ""), outer = TRUE)
        }else{
            title(paste("Stichprobenumfang: ", Vplan[i, 4], "/", Vplan[i, 5], " (Gruppe 0/Gruppe 1),", 
                "\n Zensierungen: ", Vplan.Zens[i], ", Bindungen: ", Vplan.Bdg[i], sep = ""), outer = TRUE)
        }
        dev.off()
    }
    
    
    
    ### Zeichnung der Grafiken für je zwei bis drei der Szenarien in einer getrennten Datei ###########################
    # Öffnen eines pdfs
    pdf(paste(file.erg2, "/SimEstHR_KIs_Alle.pdf", sep = ""), width = 21/2.54, height = 29.7/2.54)
      
    ## Durchlauf der einzelnen Abbildungen
    for(i in sz.reihf) {
        
        ## allgemeine Einstellungen für das Layout der Abbildungen
        layout(matrix(c(3, 3, 1:2, 6, 6, 4:5, 9, 9, 7:8), ncol = 2, byrow = TRUE), heights = c(0.1, 1, 0.1, 1, 0.1, 1))
        par(lwd = 1.2, cex = 0.6, cex.axis = 1.4, cex.lab = 1.4, cex.main = 1.45, 
            mar = c(4.5, 5, 0.5, 1), oma = c(0.1, 0.1, 3.5, 0.1))   
    
        ## Erstellung der Szenarioauswahl für eine Abbildung sowie
        ## Erzeugen der Überschriften der zwei bis drei Grafiken in einer Abbildung
        if(i < 49) {
            sz.ausw <- i + c(0, 16, 32)
            main.text <- c("Exp(variierend)/Exp(0.0693) (Gruppe 0/Gruppe 1)", 
                "Weib(0.5, variierend)/Weib(0.5, 0.2192) (Gruppe 0/Gruppe 1)", 
                "Gomp(0.3, variierend)/Gomp(0.3, 0.0109) (Gruppe 0/Gruppe 1)")
        }else{ 
            sz.ausw <- c(49, 50)
            main.text <- c("Weib(0.5, variierend)/Weib(0.5, 0.2192) (Gruppe 0/Gruppe 1), Störeffekt: 2", 
                "Weib(0.5, variierend)/Weib(0.5, 0.2192) (Gruppe 0/Gruppe 1), Störeffekt: 4")
        }
        
        ## Setzen eines Zählers für die betrachteten Szenarien innerhalb einer Abbildung 
        zaehler <- 0    
        
        ## Definition einer Liste, die alle zu plottenden Werte einer Abbildung tragen soll
        L <- vector(mode = "list", length = length(sz.ausw))
        
        ## Durchlauf der zwei bzw. drei Szenarien einer Szenarioauswahl für eine Abbildung
        for(j in sz.ausw) {
            
            # Laden der Datei des entsprechenden Szenarios
            load(paste(file.erg1, "/SimEstHR_Sz", idx.sz[j], ".RData", sep = ""))
            
            # Extrahierung der Konfidenzintervallgrenzen der drei Intervalltypten aller Hazard-Ratios eines Szenarios
            # und Berechnung der zugehörigen Überdeckungsraten und durchschnittlichen Breiten der Intervalle
            zaehler <- zaehler + 1 
            for(l in seq(along = haz.rat)) {
                konf <- sim.erg[[l]][, 14:19]
        
                cox.ud[l] <- mean((konf$Cox.KI.u < log.haz.rat[l]) & (konf$Cox.KI.o > log.haz.rat[l]), na.rm = TRUE)
                was.ud[l] <- mean((konf$Wassmer.KI.u < log.haz.rat[l]) & (konf$Wassmer.KI.o > log.haz.rat[l]), 
                    na.rm = TRUE)
                sl.ud[l] <- mean((konf$SL.KI.u < log.haz.rat[l]) & (konf$SL.KI.o > log.haz.rat[l]), na.rm = TRUE)
        
                cox.b[l] <- mean(konf$Cox.KI.o - konf$Cox.KI.u, na.rm = TRUE)
                was.b[l] <- mean(konf$Wassmer.KI.o - konf$Wassmer.KI.u, na.rm = TRUE)
                sl.b[l] <- mean(konf$SL.KI.o - konf$SL.KI.u, na.rm = TRUE)
            }
            L[[zaehler]] <- list(cox.ud = cox.ud, was.ud = was.ud, sl.ud = sl.ud, 
                cox.b = cox.b, was.b = was.b, sl.b = sl.b)
        }
        
        ## Zeichnen der Abbildungen mit den 4-6 Grafiken
        zaehler <- 0        
        for(j in sz.ausw) {
            zaehler <- zaehler + 1  
            
            # Extrahierung der bereits berechneten zu plottenden Werte eines Szenarios aus der Szenarioauswahl
            cox.ud <- L[[zaehler]]$cox.ud
            was.ud <- L[[zaehler]]$was.ud
            sl.ud <- L[[zaehler]]$sl.ud
            cox.b <- L[[zaehler]]$cox.b
            was.b <- L[[zaehler]]$was.b
            sl.b <- L[[zaehler]]$sl.b
            
            # Berechnung des zu plottenden Bereichs der y-Achse in der aktuellen Grafik der Überdeckungsraten
            min.ud <- max(0.85, min(c(0.915, sapply(L, "[[", "cox.ud"), sapply(L, "[[", "was.ud"), 
                sapply(L, "[[", "sl.ud")), na.rm = TRUE))
            max.ud <- 0.985
            
            # Zeichnen der Grafik der Überdeckungsraten für das aktuelle Szenario
            plot(log.haz.rat, cox.ud, type = "b", ylim = c(min.ud, max.ud), xlim = c(-1.25, 1.25), xaxt = "n",
                xlab = expression(Wahres ~~ ln(omega)), ylab = "Überdeckungsrate", main = "", lwd = 1.5)
            axis(1, at = log.haz.rat, labels = round(log.haz.rat, 2))
            
            lines(log.haz.rat, was.ud, type = "b", lty = "dashed", lwd = 1.5)
            lines(log.haz.rat, sl.ud, type = "b", lty = "dotted", lwd = 1.8)
            
            # Einzeichnen einer Hilfslinie für das Konfidenzniveau und einer für die Überdeckungserfolgsgrenze
            abline(h = 0.95)
            abline(h = ud.grenze, lty = "dotted")   
            
            # Erstellen einer Legende für die Niveaugrenze und die Überdeckungserfolgsgrenze
            if((j <= 16) || j == 49) {
                legend("topleft", horiz = TRUE, legend = c("Überdeckungserfolgsgrenze .9464",
                    expression("Niveau" ~ alpha == ".95")), lty = c("dotted", "solid"), lwd = 1.2, cex = 1.1)
            } 
            
            # Punkte, die außerhalb des zu plottenden y-Bereiches liegen,
            # werden als entsprechender Text in den Grafiken eingetragen.
            min.g <- par("usr")[3]
            if(!(j %in% seq(2, 48, by = 2))) {
                if(any(cox.ud < min.g))
                    text(log.haz.rat[cox.ud < min.g], 0.900, 
                        as.expression(lapply(round(cox.ud[cox.ud < min.g], 2), function(x)
                            substitute(KI[C] * ": " * s, list(s = x)))), 
                        cex = 1.1)
                if(any(was.ud < min.g))
                    text(log.haz.rat[was.ud < min.g], 0.875, 
                        as.expression(lapply(round(was.ud[was.ud < min.g], 2), function(x)
                            substitute(KI[W] * ": " * s, list(s = x)))), 
                        cex = 1.1)
                if(any(sl.ud < min.g))
                    text(log.haz.rat[sl.ud < min.g], 0.850,
                        as.expression(lapply(round(sl.ud[sl.ud < min.g], 2), function(x)
                            substitute(KI[L] * ": " * s, list(s = x)))), 
                        cex = 1.1)
            }
            else{
                if(any(cox.ud < min.g))
                    text(log.haz.rat[cox.ud < min.g], 0.900, 
                        as.expression(lapply(round(cox.ud[cox.ud < min.g], 2), function(x)
                            substitute(KI[C]^B * ": " * s, list(s = x)))), 
                        cex = 1.1)
                if(any(was.ud < min.g))
                    text(log.haz.rat[was.ud < min.g], 0.875, 
                        as.expression(lapply(round(was.ud[was.ud < min.g], 2), function(x)
                            substitute(KI[W]^B * ": " * s, list(s = x)))), 
                        cex = 1.1)
                if(any(sl.ud < min.g))
                    text(log.haz.rat[sl.ud < min.g], 0.850,
                        as.expression(lapply(round(sl.ud[sl.ud < min.g], 2), function(x)
                            substitute(KI[L]^B * ": " * s, list(s = x)))), 
                        cex = 1.1)
            }
            
            # Berechnung des zu plottenden Bereichs der y-Achse in der aktuellen Grafik der Konfidenzintervallbreiten
            min.b <- min(c(sapply(L, "[[", "cox.b"), sapply(L, "[[", "was.b"), sapply(L, "[[", "sl.b")), na.rm = TRUE)
            max.b <- max(c(sapply(L, "[[", "cox.b"), sapply(L, "[[", "was.b"), sapply(L, "[[", "sl.b")), na.rm = TRUE)
            min.b <- min.b - 0.1 * (max.b - min.b)
            
            # Zeichnen der Grafik der durchschnittlichen Konfidenzintervallbreiten für das aktuelle Szenario
            plot(log.haz.rat, cox.b, type = "b", ylim = c(min.b, max.b), xaxt = "n",
                xlab = expression(Wahres ~~ ln(omega)), ylab = "Durchschnittliche Breite", main = "", lwd = 1.5)
            axis(1, at = log.haz.rat, labels = round(log.haz.rat, 2))
            
            lines(log.haz.rat, was.b, type = "b", lty = "dashed", lwd = 1.5)
            lines(log.haz.rat, sl.b, type = "b", lty = "dotted", lwd = 1.8)
            
            # Erstellen einer Legende für die drei abbgebildeten Intervalltypen
            if((j > 32 && j < 49) || j == 50) {
                legend("bottom", xjust = 0.5, yjust = 0, horiz = TRUE, legend = c("Cox", "Wassmer", "SL"), 
                    lty = c("solid", "dashed", "dotted"), lwd = 1.8, cex = 1.2)
            }
            
            # Beschriftung der beiden Grafiken für das aktuelle Szenario    
            op <- par(mar = c(0, 0, 2.5, 0), xpd = NA)
            plot.new()
            title(main = main.text[zaehler])
            par(op)
        }
         
        ## Erstellen eines Titels für die aktuelle Abbildung  
        if(i < 49) {
            title(paste("Stichprobenumfang: ", Vplan[i, 4], "/", Vplan[i, 5], " (Gruppe 0/Gruppe 1),", 
                "\n Zensierungen: ", Vplan.Zens[i], ", Bindungen: ", Vplan.Bdg[i], ", Störeffekt: ", 
                Vplan.Stoer[i], sep = ""), outer = TRUE)
        }else{
            title(paste("Stichprobenumfang: ", Vplan[i, 4], "/", Vplan[i, 5], " (Gruppe 0/Gruppe 1),", 
                "\n Zensierungen: ", Vplan.Zens[i], ", Bindungen: ", Vplan.Bdg[i], sep = ""), outer = TRUE)
        }
    }
    dev.off()
}