# 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 Durchführung der Simulationsstudie mit allen Szenarien (ungleiche Zensierungsverteilungen) ####################################################################################################################### Simstudie.unglZens <- function(file.erg) { # Setzen eines Seeds für reproduzierbare Ergebnisse set.seed(14032012) # Spezifizierung der Anzahl der Durchläufe für jedes Szenario runs <- 10000 # Spezifizierung der Anzahl der zu betrachtenden Szenarien anz.sz <- 50 # Spezifizierung der zu betrachtenden Hazard-Ratios in jedem Szenario haz.rat <- c(1/3, 1/2, 2/3, 1/1.2, 1, 1.2, 1.5, 2, 3) ### Spezifizierung des Versuchsplans für die ersten 48 Szenarien ## Matrix der zu betrachtenden Verteilungen Dist <- data.frame(Verteilung = c("Exp", "Weib", "Gomp"), Formparameter = c(1, 0.5, 0.3), Skalenparameter = c(.06931444, .2192088, .01088074), stringsAsFactors = FALSE) ## Matrix der zu betrachtenden Stichprobengrößen Samp.Size <- matrix(c(20, 20, 100, 100, 500, 500, 50, 100), byrow = TRUE, nrow = 4) ## Vektor der zu betrachtenden Zensierungsanteile Censoring <- c(FALSE, .4) ## Rundungsgrade für die Erzeugung von Bindungen Ties <- matrix(rep(c(NA, .5), c(3, 3)), nrow = 3) ## Zur Übersicht Kombinierung der Angaben zum Versuchplan in einem Datensatz anz.sz.v <- numeric(anz.sz) Vplan <- data.frame(Verteilung = anz.sz.v, Shape = anz.sz.v, Scale = anz.sz.v, Stipro.G0 = anz.sz.v, Stipro.G1 = anz.sz.v, Zens = anz.sz.v, Bdg = anz.sz.v, Stoergr.Effekt = anz.sz.v) ### Definition eines Datensatzes zur Speicherung von Fehlermeldungen bei der Berechnung der Schätzer Error <- data.frame(Sz = 0, Haz.Rat = 0, Run = 0, Cox.S = 0, Breslow.S = 0, Wassmer.S = 0, WBC1.S = 0, WBC2.S = 0, SL.S = 0) ### Definition eines Datensatzes zur Speicherung der Szenarien, Hazard-Ratios und Durchläufe, ### in denen ein Schätzer nicht berechnet werden konnte sim.erg <- list() ### Definition eines Zählers für die Szenarien ### und einen für die Durchläufe mit Fehlermeldung in mindestens einem Schätzer z1 <- 1; z2 <- 1 idx.run <- 1:runs ### Durchführung der Simulationsstudie ## Durchlauf aller spezifizierten Verteilungen for(i in 1:dim(Dist)[1]) { distrib <- Dist[i, 1] shape <- rep(Dist[i, 2], 2) Scal <- matrix(c(haz.rat * Dist[i, 3], rep(Dist[i, 3], 9)), ncol = 2) ## Durchlauf aller spezifizierten Stichprobengrößen for(j in 1:dim(Samp.Size)[1]) { samp.size <- Samp.Size[j, ] n0 <- samp.size[1] n1 <- samp.size[2] n <- n0 + n1 ## Durchlauf aller spezifizierten Zensierungsanteile for(k in seq(along = Censoring)) { censoring <- Censoring[k] ## Berechnung des Skalenparameters der Exponentialverteilung der Zensierungsverteilung für Gruppe 1, ## der für jedes Hazard-Ratio gleich ist, wenn Zensierungen vorgesehen sind ## (vgl. Datei 'Sim_OptZens.R') if(censoring) { cens.par.1 <- optimize(opt.cens.par.1, interval = c(0, 500), cens.rate = censoring, distrib = distrib, shape = shape[1], scal = Dist[i, 3])$minimum } else{ cens.par.1 <- 0 } ## Durchlauf aller spezifizierten Rundungsgrade for(l in 1:dim(Ties)[2]) { ties <- Ties[i, l] ## Durchlauf aller spezifizierten Hazard-Ratios for(m in seq(along = haz.rat)){ scal <- Scal[m, ] ## Berechnung des Skalenparameters der Exponentialverteilung der Zensierungsverteilung ## für Gruppe 2, der für jedes Hazard-Ratio anders ist, wenn Zensierungen vorgesehen sind ## (vgl. Datei 'Sim_OptZens.R') if(censoring) { cens.par.0 <- optimize(opt.cens.par.1, interval = c(0.000001, 100), cens.rate = censoring, distrib = distrib, shape = shape[1], scal = scal[1])$minimum } else{ cens.par.0 <- 0 } ## Simulation für das spezifizierte Szenario und Hazard-Ratio sim.erg.temp <- Sim(runs = runs, distrib = distrib, shape = shape, scal = scal, samp.size = samp.size, censoring = censoring, cens.par.0 = cens.par.0, cens.par.1 = cens.par.1,ties = ties, stoergr = FALSE, stoergr.corr = FALSE) ## Übertragung der Ergebnisse eines Szenarios und Hazard-Ratios in die oben definierte Liste sim.erg[[m]] <- sim.erg.temp ## Speicherung der Szenarien, Hazard-Ratios, Durchläufe und Schätzer, ## für die kein Schätzwert berechnet werden konnte sim.erg.est <- sim.erg.temp[, 1:6] idx.error <- apply(sim.erg.est, 1, function(x) any(is.na(x))) if(any(idx.error)) { idx.run.error <- idx.run[idx.error] for(o in idx.run.error) { anzahl <- is.na(sim.erg.est[o, ]) Error[z2, ] <- c(z1, m, o, anzahl) z2 <- z2 + 1 } } } ## Speicherung der Eigenschaften der Szenarien im Übersichtsdatensatz 'Vplan' Vplan[z1, ] <- c(distrib, shape[1], scal[2], samp.size, ifelse(censoring, "40%", "keine"), ifelse(is.na(ties), "keine", "siehe Grafik"), "keiner") ## Speicherung der Ergebnisse für ein Szenario (alle Durchläufe und Hazard-Ratios) ## in einer 'RData'-Datei save(sim.erg, file = paste(file.erg, "/SimEstHR_Sz", z1, ".RData", sep = "")) z1 <- z1 + 1 } } } } ### Spezifizierung und Durchführung der zwei zusätzlichen Szenarien mit einer Störgröße ## Definition der Eigenschaften i <- 2 distrib <- Dist[i, 1] shape <- rep(Dist[i, 2], 2) Scal <- matrix(c(haz.rat * Dist[i, 3], rep(Dist[i, 3], 9)), ncol = 2) samp.size <- c(100, 100) censoring <- FALSE ties <- NA stoergr <- c(2, 4) ## Durchführung der zwei Szenarien for(p in 1:2) { for(m in seq(along = haz.rat)){ scal <- Scal[m, ] sim.erg.temp <- Sim(runs = runs, distrib = distrib, shape = shape, scal = scal, samp.size = samp.size, censoring = censoring, ties = ties, stoergr = stoergr[p], stoergr.corr = TRUE) sim.erg[[m]] <- sim.erg.temp sim.erg.est <- sim.erg.temp[, 1:6] idx.error <- apply(sim.erg.est, 1, function(x) any(is.na(x))) if(any(idx.error)) { idx.run.error <- idx.run[idx.error] for(o in idx.run.error) { anzahl <- is.na(temp.na.est[o, ]) Error[z2, ] <- c(z1, m, o, anzahl) z2 <- z2 + 1 } } } ## Speicherung der Eigenschaften der zwei Szenarien im Übersichtsdatensatz 'Vplan' Vplan[z1, ] <- c(distrib, shape[1], scal[2], samp.size, "keine", "keine", stoergr[p]) ## Speicherung der Ergebnisse der zwei Szenarien in einer 'R-Data'-Datei save(sim.erg, file = paste(file.erg, "/SimEstHR_Sz", z1, ".RData", sep = "")) z1 <- z1 + 1 } ### Speicherung des Datensatzes über den Versuchsplan und den der Fehlerstatistik zu den einzelnen Schätzern save(Vplan, file = paste(file.erg, "/SimEstHR_Info_Vplan.RData", sep = "")) save(Error, file = paste(file.erg, "/SimEstHR_Info_Error.RData", sep = "")) }