#  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 Berechnung der Punkt- und Intervallschätzer für alle Durchläufe eines gegebenen Szenarios
#######################################################################################################################
Sim <- function(runs = 100, distrib = c("Exp", "Weib", "Gomp"), shape = c(1, 1), scal = c(1, 1), 
    samp.size = c(20, 20), censoring = FALSE, cens.par.0 = 0.4, cens.par.1 = 0.4, ties = NA, 
    stoergr = FALSE, stoergr.corr = FALSE) { 
    
    distrib <- match.arg(distrib)
    
    
    ### Definition der Vektoren, die von der Funktion ausgegeben werden
    ### und für jeden Simulationsdurchlauf eine Information tragen
    ## einseitige Teststatistik des Logrank-Tests
    LR.Stat.v <- numeric(runs)
    LR.Stat.v[] <- NA
    
    ## Schätzwerte der sechs zu untersuchenden Punktschätzer
    Cox.S <- Breslow.S <- Wassmer.S <- WBC1.S <- WBC2.S <- SL.S <- LR.Stat.v
    
    ## Anzahl an Bindungen unter den unzensierten Überlebenszeiten sowie Anzahl der Ereignisse 
    ## jeweils gruppenübergreifend sowie für Gruppe 0 und 1 getrennt in der zu ziehenden Stichprobe
    d.v <- d0.v <- d1.v <- bd0.v <- bd1.v <- bd01.v <- LR.Stat.v 
    
    ## Untere/obere Intervallgrenzen der drei zu untersuchenden Intervallschätzer
    Cox.KI.u <- Cox.KI.o <- Wassmer.KI.u <- Wassmer.KI.o <- SL.KI.u <- SL.KI.o <-LR.Stat.v
        
        
    ### Extrahierung des Stichprobenumfangs und -verhältnisses
    n0 <- samp.size[1] 
    n1 <- samp.size[2]
    n <- n0 + n1
    r <- n0 / n1
        
    
    ### Durchführung der Simulationen für das gegebene Szenario 
    ### (ein Schleifendurchlauf entspricht einem Simulationsdurchlauf)
    for(i in 1:runs) {
        
        ### Erzeugung des Effekts und der Identität der Störgröße für alle n Patienten, wenn spezifiziert
        if(!stoergr) { stoergr.eff <- rep(1, n)
        } else { 
            stoergr.id <- rbinom(n, 1, 1/2)
            stoergr.eff <- ifelse(stoergr.id == 0, 1, stoergr)
        }
        
        
        ### Erzeugung der Überlebenszeiten für alle n Patienten
        if((distrib == "Exp") || (distrib == "Weib")) {
            etime <- rweibull(n, shape = rep(shape, samp.size), 
                scale = (rep(scal, samp.size) * stoergr.eff)^(-1/rep(shape, samp.size)))
        } else {
            etime <- rgompertz(n, alpha = rep(shape, samp.size), theta = rep(scal, samp.size) * stoergr.eff)        
        } 
        
        
        ### Notierung der Gruppenzugehörigkeit aller n Patienten (1 und 2 entspricht 0 bzw. 1 in der Arbeit)
        group <- rep(1:2, samp.size)
        

        ### Erzeugung der Zensierungszeiten, Notierung des Zensierungsstatus 
        ### und entsprechende Reduktion der (unzensierten) Überlebenszeiten aller n Patienten
        if(!censoring) {
            status <- rep(1, n) 
        } else {
            z0 <- rexp(n0, cens.par.0)
            z1 <- rexp(n1, cens.par.1)
            z <- c(z0, z1)
            status <- as.numeric(etime <= z)
            etime <- ifelse(etime <= z, etime, z)
        }
        
                    
        ### Erzeugung von Bindungen, wenn spezifiziert (wie in Abschnitt 4.2.1 beschrieben)
        if(!is.na(ties)) {
            if(ties >= 1) {
                etime <- floor(etime * 10^ties) / 10^ties
            }
            else {
                etime <- etime %/% ties / (1 / ties)
            } 
        }
        
        
        ### Bestimmung der Anzahl der gebundenen unzensierten Überlebenszeiten
        ### nach Gruppe getrennt sowie gruppenübergreifend
        idx.n0 <- 1:n0    
        bd0 <- sum(duplicated(etime[idx.n0][status[idx.n0] == 1])) 
        bd0.v[i] <- bd0
        bd1 <- sum(duplicated(etime[-idx.n0][status[-idx.n0] == 1]))
        bd1.v[i] <- bd1 
        bd01 <- sum(duplicated(etime[status == 1]))
        bd01.v[i] <- bd01
        
        
        ### Berechnung der Schätzwerte der sechs Punktschätzer
        ## Cox-Schätzer
        group.cox <- rep(1:0, samp.size)
        opw <- options(warn = 2)
        
        # Im Cox-Modell kann die Störgröße mitmodelliert werden
        if(!stoergr.corr) { 
            cox.mo <- try(coxph(Surv(etime, status) ~ group.cox))
        } else { 
            cox.mo <- try(coxph(Surv(etime, status) ~ group.cox + stoergr.id))
        }
        
        # In extremen Szenarien kann es passieren, dass kein Schätzwert für den Cox-Schätzer bestimmt werden kann
        if(inherits(cox.mo, "try-error")){
            Cox.S.temp <- NA
        } else {
            Cox.S.temp <- as.numeric(cox.mo$coef[1])
        }
        Cox.S[i] <- Cox.S.temp
        
        
        ## Berechnung der Logrank-Statistik, die für die Berechnung des Wassmer, WBC-I-, WBC-II- und SL-Schätzers 
        ## benötigt wird. Der Code der beiden Funktionen 'LR.oB()' und LR() für den Fall ohne und mit Bindungen 
        ## ist in der Datei 'Sim_LRStatistiken.R' zu finden.
        if(is.na(ties) & all(!duplicated(etime))) {
            LR.temp <- try(LR.oB(etime = etime, group = group, status = status))
        } else {
            LR.temp <- try(LR(etime = etime, group = group, status = status))
        }   
        options(opw)   
        
        # Analog zur Berechnung des Cox-Modells sollen die Simulationen weiterlaufen, 
        # wenn in einem Durchlauf bei der Berechnung der Logrank-Statistik eine Fehlermeldung auftritt.
        if(inherits(LR.temp, "try-error")) next
        if(is.na(LR.temp$LR.T)) next
        
        # Extrahierung der einseitigen Teststatistik des Logrank-Tests
        LR.Stat <- LR.temp$LR.T; LR.Stat.v[i] <- LR.Stat
        
        # Extrahierung der Anzahl der Ereignisse (gruppenübergreifend sowie für Gruppe 0 und 1 getrennt)
        d <- LR.temp$d
        d.v[i] <- d
        d0 <- LR.temp$d0
        d0.v[i] <- d0
        d1 <- (d - d0)
        d1.v[i] <- d1
        
        # Extrahierung der Anzahl der Ereignisse, Anzahl der Patienten unter Risiko 
        # und erwarteter Anzahl der Ereignisse an den einzelnen Ereigniszeiten einer Stichprobe
        d_1j <- LR.temp$d_1j
        d_2j <- LR.temp$d_2j
        n_1j <- LR.temp$n_1j
        n_2j <- LR.temp$n_2j
        e_1j <- LR.temp$e_1j
        e_2j <- LR.temp$e_2j
        
        # Bestimmung der Anzahl der Ereigniszeiten
        if(!is.na(ties)) {delta <- LR.temp$delta}
        
        
        ## Breslow-Schätzer
        # Beim Breslow-Schätzer sollen die Fälle 'abgefangen' werden, in denen keine Ereignisse 
        # in einer Gruppe auftreten, womit dann die Werte 'Inf' bzw. '-Inf' als Schätzwerte resultieren.
        opw <- options(warn = 2)
        Breslow.S.temp <- try(log((d0 * sum(e_2j))  / (d1 * sum(e_1j))))
        if(!inherits(Breslow.S.temp, "try-error") & (Breslow.S.temp != Inf) & 
            (Breslow.S.temp != -Inf)) {
            Breslow.S[i] <- Breslow.S.temp
        }
        options(opw)  
        
        
        ## Wassmer-Schätzer
        r.est <- (1 + r) / sqrt(r)
        Wassmer.S.temp <- r.est / sqrt(d) * LR.Stat
        Wassmer.S[i] <- Wassmer.S.temp
                
        
        ## Bias-Corrected-Wassmer-Schätzer I
        var.LR <- LR.temp$var.LR
        bias1 <- r.est / sqrt(d) / sqrt(var.LR) * 
            sum(n_1j * exp(Wassmer.S.temp) * (d_1j + d_2j) / (n_1j * 
                exp(Wassmer.S.temp) + n_2j) - e_1j) - Wassmer.S.temp
        WBC1.S.temp <- Wassmer.S.temp - bias1
        WBC1.S[i] <- WBC1.S.temp
        
        
        ## Bias-Corrected-Wassmer-Schätzer II
        n_1j.exp <- n0
        n_2j.exp <- n1
        if(LR.Stat < 0) {lhr <- log(1/2)} else {lhr <- log(2)}

        if(is.na(ties)) {
            for(j in 1:(d - 1)) {
                p_1j <- n_1j.exp[j] * exp(lhr) / (n_1j.exp[j] * exp(lhr) + n_2j.exp[j])
                n_1j.exp[j + 1] <- n_1j.exp[j] - p_1j
                n_2j.exp[j + 1] <- n_2j.exp[j] - (1 - p_1j)
            }
            n_j.exp <- n_1j.exp + n_2j.exp
            var.LR.exp <- sum((n_1j.exp * n_2j.exp) / n_j.exp^2)
        
            bias2 <- r.est / sqrt(d) / sqrt(var.LR.exp) * 
                sum(n_1j.exp * exp(lhr) / (n_1j.exp * exp(lhr) + n_2j.exp) - n_1j.exp / n_j.exp) - lhr
        } else {        
            d.est <- d / delta
            for(j in 1:(delta - 1)) {
                p_1j <- n_1j.exp[j] * exp(lhr) / (n_1j.exp[j] * exp(lhr) + n_2j.exp[j])
                n_1j.exp[j + 1] <- n_1j.exp[j] - p_1j * d.est
                n_2j.exp[j + 1] <- n_2j.exp[j] - (1 - p_1j) * d.est
            }
            n_j.exp <- n_1j.exp + n_2j.exp
            d.est.delta <- rep(d.est, delta)
            var.LR.exp <- sum((n_1j.exp * n_2j.exp * d.est.delta * (n_j.exp - d.est.delta)) / 
                (n_j.exp^2 * (n_j.exp - 1)))
        
            bias2 <- r.est / sqrt(d) / sqrt(var.LR.exp) * 
                sum(n_1j.exp * exp(lhr) * d.est / (n_1j.exp * exp(lhr) + n_2j.exp) - 
                n_1j.exp * d.est / n_j.exp) - lhr
        }
        WBC2.S.temp <- Wassmer.S.temp - bias2
        WBC2.S[i] <- WBC2.S.temp
                
        
        ## SL-Schätzer
        SL.S.temp <- r.est^2 / d * (d0 - sum(e_1j))
        SL.S[i] <- SL.S.temp
        
        
        ### Berechnung der unteren u. oberen Grenzen der drei Konfidenzintervalle
        ## Cox-Konfidenzintervall
        if(inherits(cox.mo, "try-error")){
            Cox.KI.u[i] <- NA
            Cox.KI.o[i] <- NA
        } else {
            Cox.KI.u[i] <- Cox.S.temp - sqrt(cox.mo$var[1,1]) * qnorm(0.975)
            Cox.KI.o[i] <- Cox.S.temp + sqrt(cox.mo$var[1,1]) * qnorm(0.975)
        }
        
        
        ## Wassmer-Konfidenzintervall
        p_1j.w <- n_1j * exp(Wassmer.S.temp) / 
            (n_1j * exp(Wassmer.S.temp) + n_2j)
        p_2j.w <- n_2j / (n_1j * exp(Wassmer.S.temp) + n_2j)
        
        Wassmer.KI.u[i] <- Wassmer.S.temp - r.est / sqrt(d) * 
            sqrt(sum((d_1j + d_2j) * p_1j.w * p_2j.w)/ var.LR) * qnorm(0.975)
        Wassmer.KI.o[i] <- Wassmer.S.temp + r.est / sqrt(d) * 
            sqrt(sum((d_1j + d_2j) * p_1j.w * p_2j.w)/ var.LR) * qnorm(0.975)
        
        ## SL-Konfidenzintervall
        p_1j.v <- n_1j * exp(SL.S.temp) / (n_1j * exp(SL.S.temp) + n_2j)
        p_2j.v <- n_2j / (n_1j * exp(SL.S.temp) + n_2j)
        
        SL.KI.u[i] <- SL.S.temp - (r.est)^2 / d * sqrt(sum((d_1j + d_2j) * 
            p_1j.v * p_2j.v)) * qnorm(0.975)
        SL.KI.o[i] <- SL.S.temp + (r.est)^2 / d * sqrt(sum((d_1j + d_2j) * 
            p_1j.v * p_2j.v)) * qnorm(0.975)
    }
    
    
    ### Ausgabe der Ergebnisse
    return(data.frame(Cox.S = Cox.S, Breslow.S = Breslow.S, 
        Wassmer.S = Wassmer.S, WBC1.S = WBC1.S, WBC2.S = WBC2.S, SL.S = SL.S, 
        LR.Stat.v = LR.Stat.v, d.v = d.v, d0.v = d0.v, d1.v = d1.v, bd0.v = bd0.v, bd1.v = bd1.v, bd01.v, 
        Cox.KI.u = Cox.KI.u, Cox.KI.o = Cox.KI.o, 
        Wassmer.KI.u = Wassmer.KI.u, Wassmer.KI.o = Wassmer.KI.o, SL.KI.u = SL.KI.u, SL.KI.o = SL.KI.o))
}