#  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.


#######################################################################################################################
### Funktionen zur Berechung der einseitigen Teststatistik des Logrank-Tests
#######################################################################################################################

###################################################################################################
### a) Fall ohne Bindungen
###################################################################################################
LR.oB <- function(etime, group, status) {
                            
    ### Erzeugen der Anzahl Objekte unter Riskio pro Ausfallzeitpunkt
    tab.ti.gr <- table(etime, group)
    l.etime <- length(etime)
    
    n1 <- sum(group == 1)
    n2 <- sum(group == 2)
    n_1j <- n1 - c(0, cumsum(as.numeric(tab.ti.gr[-l.etime, 1])))
    n_2j <- n2 - c(0, cumsum(as.numeric(tab.ti.gr[-l.etime, 2])))
        
    if(any(status == 0)) {
        tab.ti.st <- table(etime, status)
        index.nc <- tab.ti.st[, 2] > 0
        n_1j <- n_1j[index.nc]
        n_2j <- n_2j[index.nc]
    }
    
    n_j <- n_1j + n_2j                                                                   
                 
                            
    ### Erzeugen der Anzahl Ausfälle pro Ausfallzeitpunkt
    group.nc <- group[status == 1]
    etime.nc <- etime[status == 1]
    tab.ti.gr.nc <- table(etime.nc, group.nc)  
    
    if(any(group.nc == 1))            
        d_1j <- as.numeric(tab.ti.gr.nc[, "1"])
    else 
        d_1j <- rep(0, length(n_j))
        
    d_2j <- 1 - d_1j 
    d_j <- d_1j + d_2j
    d <- sum(status == 1)
    d0 <- sum(d_1j)    
    
    
    ### Erzeugen der Anzahl erwarteter Ausfälle pro Ausfallzeitpunkt
    e_1j <- n_1j / n_j
    e_2j <- n_2j / n_j
    
    
    ### Berechnung der Teststatistik
    var.LR <- sum((n_1j * n_2j) / n_j^2)
    LR.T <- sum(d_1j - e_1j) / sqrt(var.LR)
    
    
    ### Ausgabe der Teststatistik und anderer berechneter Größen
    return(list(LR.T = LR.T, d = d, d0 = d0, d_1j = d_1j, d_2j = d_2j, d_j = d_j,
        n_1j = n_1j, n_2j = n_2j, e_1j = e_1j, e_2j = e_2j, var.LR = var.LR))
}



###################################################################################################
### b) allgemeiner Fall (mit und ohne Bindungen)
###################################################################################################
LR <- function(etime, group, status) {
    
    ### Vorbereitung: Wenn der größte Eintrag in 'etime' eindeutig und nicht-zensiert ist, 
    ### ergibt sich für die größte Ereigniszeit 'n_{.j} = 1' und damit "(n_{.j} - 1) = 0".
    ### Durch letzteren Term wird in der Varianzformel in einem Summanden dividiert.
    ### Damit dieses Problem nicht auftritt, wird die größte Beobachtung zensiert, 
    ### wenn sie eindeutig ist.
    etime.nc <- etime[status == 1] 
    max.etime.nc <- max(etime.nc)  

    if((max.etime.nc == max(etime)) && (sum(max.etime.nc == etime) == 1)) {
        status[which.max(etime)] <- 0
    }  
        
                        
    ### Erzeugen der Anzahl Objekte unter Riskio pro Ausfallzeitpunkt
    tab.ti.gr <- table(etime, group)
    l.etime.u <- length(unique(etime))
    
    n1 <- sum(group == 1)
    n2 <- sum(group == 2)
    n_1j <- n1 - c(0, cumsum(as.numeric(tab.ti.gr[-l.etime.u, 1])))
    n_2j <- n2 - c(0, cumsum(as.numeric(tab.ti.gr[-l.etime.u, 2])))
        
    if(any(status == 0)) {
        tab.ti.st <- table(etime, status)
        index.nc <- tab.ti.st[, 2] > 0
        n_1j <- n_1j[index.nc]
        n_2j <- n_2j[index.nc]
    }
    
    n_j <- n_1j + n_2j                                   
                     
                            
    ### Erzeugen der Anzahl Aufälle pro Ausfallzeitpunkt 
    group.nc <- group[status == 1]
    etime.nc <- etime[status == 1]
    tab.ti.gr.nc <- table(etime.nc, group.nc)  
    
    if(any(group.nc == 1))            
        d_1j <- as.numeric(tab.ti.gr.nc[, "1"])
    else 
        d_1j <- rep(0, length(n_j))
        
    if(any(group.nc == 2))            
        d_2j <- as.numeric(tab.ti.gr.nc[, "2"])
    else 
        d_2j <- rep(0, length(n_j))
    
    d_j <- d_1j + d_2j 
    d0 <- sum(d_1j)
    
    
    ### Erzeugen der Anzahl erwarteter Ausfälle pro Ausfallzeitpunkt
    e_1j <- n_1j * d_j / n_j
    e_2j <- n_2j * d_j / n_j
    
    
    ### Berechnung der Teststatistik
    var.LR <- sum((n_1j * n_2j * d_j * (n_j - d_j)) / ((n_j^2) * (n_j - 1)))
    LR.T <- sum(d_1j - e_1j) / sqrt(var.LR)
    d <- sum(d_j)
    delta <- length(d_j)
    r <- n1 / n2
    
    
    ### Ausgabe der Teststatistik und anderer berechneter Größen
    return(list(LR.T = LR.T, d = d, d0 = d0, d_1j = d_1j, d_2j = d_2j, d_j = d_j,
        n_1j = n_1j, n_2j = n_2j, e_1j = e_1j, e_2j = e_2j, delta = delta, var.LR = var.LR))
}