# 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 parallelen Boxplots für die Ergebnisse der Punktschätzverfahren ####################################################################################################################### SimEstHR.Boxplots <- function(file.erg) { ### Vorbereitungen ## Laden des Versuchsplans load(file = paste(file.erg, "/SimEstHR_Info_Vplan.RData", sep = "")) ## Spezifikationen, die für Berechnungen oder den Titel einer Grafik benötigt werden Vplan[Vplan[, 6] == "40%", 6] <- "40% (Gruppe 0/Gruppe 1: siehe Grafik)" Samp.Size <- matrix(c(20, 20, 100, 100, 500, 500, 50, 100), byrow = TRUE, nrow = 4) 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) ## Definition von Punkten, die für das Layout der Grafiken benötigt werden fuer.at <- c(-0.25, -0.15, -0.05, 0.05, 0.15, 0.25) xaxis.2 <- unlist(lapply(seq(0, 56, by = 7), function(x) x + 1:6)) xaxis.2.re <- seq(3.5, 52.5, by = 7) max.1 <- c(3.4, 2, 1.45, 2.4) min.1 <- -max.1 min.2 <- c(-1.8, 0, 0.2, -0.2) max.2 <- c(10.2, 5, 4, 6.2) startmax1 <- c(0.3, 0.25, 0.18, 0.25) startmax2 <- c(0.7, 0.3, 0.25, 0.35) startmin1 <- 0.1 startmin2 <- 0.02 xlim1 <- c(log.haz.rat[1] + fuer.at[1] * 0.25 - 0.08, log.haz.rat[9] + fuer.at[6] * 0.25 + 0.02) xlim2 <- c(0.8, 62.1) ## Definition der Namen der Punktschätzer name <- c(expression(S[C]), expression(S[B]), expression(S[W]), expression(S[W]^{B[1]}), expression(S[W]^{B[2]}), expression(S[L])) name.B <- c(expression(S[C]^B), expression(S[B]^B), expression(S[W]^B), expression(S[W]^list(B[1],B)), expression(S[W]^list(B[2],B)), expression(S[L]^B)) name.exp <- c(expression(exp(S[C])), expression(exp(S[B])), expression(exp(S[W])), expression(exp(S[W]^{B[1]})), expression(exp(S[W]^{B[2]})), expression(exp(S[L]))) name.expB <- c(expression(exp(S[C]^B)), expression(exp(S[B]^B)), expression(exp(S[W]^B)), expression(exp(S[W]^list(B[1],B))), expression(exp(S[W]^list(B[2],B))), expression(exp(S[L]^B))) mycols <- c("blue", "green", "red", "purple", "orange", "yellow") ## Definition von Datensätzen, die für die beiden Grafiken die für ein Szenario und (log.) Hazard-Ratio ## Anzahl an Punkten außerhalb der Plotregionen (pro Schätzer) beinhalten sollen out.range.lw <- data.frame(Sz = 0, log.haz.rat = 0, Cox.S = 0, Breslow.S = 0, Wassmer.S = 0, WBC1.S = 0, WBC2.S = 0, SL.S = 0) out.range.w <- data.frame(Sz = 0, haz.rat = 0, Cox.S = 0, Breslow.S = 0, Wassmer.S = 0, WBC1.S = 0, WBC2.S = 0, SL.S = 0) yWerte.lw <- matrix(0, nrow = 50, ncol = 2) yWerte.w <- matrix(0, nrow = 50, ncol = 2) ### Zeichnung der beiden Grafiken für jedes Szenario in einer getrennten Datei #################################### ## Definition von Zählern für die Punkte, die außerhalb des jeweiligen Plotbereiches liegen ## (einmal für die Grafiken des logarithmierten Hazard-Ratios und einmal für die des Hazard-Ratios selbst) z1 <- 1 z2 <- 1 for(i in 1:50) { ## Laden der Datei des entsprechenden Szenarios und Öffnen eines pdfs load(paste(file.erg, "/SimEstHR_Sz", i, ".RData", sep = "")) pdf(paste(file.erg, "/SimEstHR_Boxplots_", i, ".pdf", sep = ""), width = 21/2.54, height = 29.7/2.54) ## erste Grafik (für das logarithmierte Hazard-Ratio) ############################################### # Extrahierung der Schätzwerte aller sechs Schätzer est <- as.data.frame(lapply(sim.erg, "[", 1:6)) # allgemeine Einstellungen für das Layout der Grafiken par(mfrow = c(2, 1), lwd = 1.1, cex = 0.6, cex.axis = 1.4, cex.lab = 1.4, cex.main = 1.4, mar = c(4.5, 5, 0.5, 1.5), oma = c(0.1, 0.1, 3.5, 0.1)) min.sz <- min.1[which(Samp.Size[, 1] == Vplan$Stipro.G0[i])] max.sz <- max.1[which(Samp.Size[, 1] == Vplan$Stipro.G0[i])] startmax.sz <- startmax1[which(Samp.Size[, 1] == Vplan$Stipro.G0[i])] # Zeichnen der parallelen Boxplots boxplot(est, at = unlist(lapply(log.haz.rat, function(x) x + fuer.at * 0.25)), xlim = xlim1, ylim = c(min.sz, max.sz), las = 1, pars = list(boxwex = 0.017, staplewex = 0.5), col = mycols, xaxt = "n", cex = 0.2, xlab = expression(Wahres ~~ ln(omega)), ylab = expression(Geschätztes ~~ ln(omega))) axis(1, at = log.haz.rat, labels = round(log.haz.rat, 2)) for(j in seq(along = haz.rat)) { lines(c(log.haz.rat[j] + fuer.at[1] * 0.25 - 0.02, log.haz.rat[j] + fuer.at[6] * 0.25 + 0.02), rep(log.haz.rat[j], 2)) } # Schreiben einer Legende für die Zuordnung der Schätzer zu den einzelnen Boxplots if(Vplan$Bdg[i] == "keine") { legend("topleft", horiz = TRUE, legend = name, col = mycols, lwd = 1.8, lty = "solid", cex = 1.2) } else { legend("topleft", horiz = TRUE, legend = name.B, col = mycols, lwd = 1.8, lty = "solid", cex = 1.2) } # Eintragen der gruppenweisen Zensierungsanteile in die Grafik, wenn Zensierungen vorliegen texth <- 0.95 sh <- strheight(expression(S[W]), cex = texth) sh.fac <- 1.5 if(Vplan$Zens[i] != "keine") { n0 <- as.numeric(Vplan$Stipro.G0[i]); n1 <- as.numeric(Vplan$Stipro.G1[i]) zens0 <- (n0 - unlist(lapply(lapply(sim.erg, "[", 9), colMeans))) / n0 * 100 zens1 <- (n1 - unlist(lapply(lapply(sim.erg, "[", 10), colMeans))) / n1 * 100 text(x = c(xlim1[1] - 0.05, log.haz.rat[5]), y = c(max.sz - startmax.sz, min.sz + startmin1) + sh * sh.fac * c(-0, 8), labels = "Zens(%): ", cex = texth, adj = 0) text(x = log.haz.rat, y = rep(c(max.sz - startmax.sz, min.sz + startmin1) + sh * sh.fac * c(-0, 8), c(5, 4)), labels = paste(formatC(zens0, digits = 0, format = "f"), "/", formatC(zens1, digits = 0, format = "f"), sep = ""), cex = texth) } # Eintragen des gruppenübergreifenden Bindungssanteils in die Grafik, wenn Bindungen vorliegen if(Vplan$Bdg[i] != "keine") { n <- as.numeric(Vplan$Stipro.G0[i]) + as.numeric(Vplan$Stipro.G1[i]) bd01 <- numeric(length(haz.rat)) for(l in seq(along = haz.rat)) { sim.erg.l <- sim.erg[[l]] bd01[l] <- mean(sim.erg.l[, 13] / sim.erg.l[, 8]) * 100 } text(x = c(xlim1[1] - 0.05, log.haz.rat[5]), y = c(max.sz - startmax.sz, min.sz + startmin1) + sh * sh.fac * c(-1, 7), labels = "Bdg(%): ", cex = texth, adj = 0) text(x = log.haz.rat, y = rep(c(max.sz - startmax.sz, min.sz + startmin1) + sh * sh.fac * c(-1, 7), c(5, 4)), labels = paste(formatC(bd01, digits = 0, format = "f"), sep = ""), cex = texth) } # Eintragen der MSEs in die Grafik für die einzelnen logarithmierten Hazard-Ratios und Schätzer sowie # Speicherung der Anzahl d. Punkte außerhalb des Plotbereichs pro Schätzer für die einzelnen log. Hazard-Ratios yWerte.lw[i, ] <- par("usr")[c(3, 4)] for(l in seq(along = haz.rat)) { est.temp <- sim.erg[[l]][,1:6] MSE <- apply(est.temp, 2, function(x) var(x, na.rm = TRUE)) + (apply(est.temp, 2, function(x) mean(x, na.rm = TRUE)) - log.haz.rat[l])^2 max.MSE <- which(MSE == max(MSE)) if(l <= 5) { text(x = log.haz.rat[l], y = max.sz - startmax.sz - sh * sh.fac * seq(3, 8, by = 1), labels = formatC(1000 * MSE, digits = 2, format = "f"), col = ifelse(MSE == min(MSE), "red", "black"), adj = 0.5, cex = texth)} else { text(x = log.haz.rat[l], y = min.sz + startmin1 + sh * sh.fac * seq(5, 0, by = -1), labels = formatC(1000 * MSE, digits = 2, format = "f"), col = ifelse(MSE == min(MSE), "red", "black"), adj = 0.5, cex = texth)} anzahl1 <- apply(est.temp, 2, function(x) sum(((x < yWerte.lw[i, 1]) | (x > yWerte.lw[i, 2])), na.rm = TRUE)) if(any(anzahl1)) { out.range.lw[z1, ] <- c(i, l, anzahl1) z1 <- z1 + 1 } } # Eintragen der MSE- und Schätzer-Beschriftungen text(x = c(xlim1[1] - 0.05, log.haz.rat[5]), y = c(max.sz - startmax.sz, min.sz + startmin1) + sh * sh.fac * c(-2, 6), labels = paste("MSE * 1000:"), adj = 0) text(x = c(rep(xlim1[1] - 0.05, 6), rep(log.haz.rat[5], 6)), y = c(max.sz - startmax.sz - sh * sh.fac * seq(3, 8, by = 1), min.sz + startmin1 + sh * sh.fac * seq(5, 0, by = -1)), labels = rep(name, 2), cex = texth, adj = 0) ## zweite Grafik (für das nicht-logarithmierte Hazard-Ratio) ######################################## # Extrahierung Schätzwerte aller sechs Schätzer estw <- exp(est) # allgemeine Einstellungen für das Layout der Grafiken par(mar = c(4.5, 5, 2.5, 1.5)) min.sz <- min.2[which(Samp.Size[, 1] == Vplan$Stipro.G0[i])] max.sz <- max.2[which(Samp.Size[, 1] == Vplan$Stipro.G0[i])] startmax.sz <- startmax2[which(Samp.Size[, 1] == Vplan$Stipro.G0[i])] # Zeichnen der parallelen Boxplots boxplot(estw, at = xaxis.2, xlim = c(0.9, 62.1), ylim = c(min.sz, max.sz), las = 1, pars = list(boxwex = 0.7, staplewex = 0.5), col = mycols, xaxt = "n", cex = 0.2, xlab = expression(Wahres ~~ omega), ylab = expression(Geschätztes ~~ omega)) axis(1, at = seq(3.5, 59.5, by = 7), labels = round(haz.rat, 2)) for(k in seq(along = haz.rat)) { lines(c(seq(1, 57, by = 7)[k] - 0.3, seq(6, 62, by = 7)[k]) + 0.3, rep(haz.rat[k], 2)) } # Schreiben einer Legende für die Zuordnung der Schätzer zu den einzelnen Boxplots if(Vplan$Bdg[i] == "keine") { legend("topleft", horiz = TRUE, legend = name.exp, col = mycols, lwd = 1.8, lty = "solid", cex = 1.2) } else { legend("topleft", horiz = TRUE, legend = name.expB, col = mycols, lwd = 1.8, lty = "solid", cex = 1.2) } # Eintragen der MSEs in die Grafik für die einzelnen Hazard-Ratios und Schätzer sowie # Speicherung der Anzahl d. Punkte außerhalb des Plotbereichs pro Schätzer für die einzelnen Hazard-Ratios sh <- strheight(expression(S[W]), cex = texth) yWerte.w[i, ] <- par("usr")[c(3, 4)] for(l in seq(along = haz.rat)) { estw.temp <- exp(sim.erg[[l]][, 1:6]) MSEw <- apply(estw.temp, 2, function(x) var(x, na.rm = TRUE)) + (apply(estw.temp, 2, function(x) mean(x, na.rm = TRUE)) - haz.rat[l])^2 if(l <= 8) { text(x = xaxis.2.re[l], y = max.sz - startmax.sz - sh * sh.fac * seq(1, 6, by = 1), labels = formatC(1000 * MSEw, digits = 2, format = "f"), col = ifelse(MSEw == min(MSEw), "red", "black"), adj = 0.5, cex = texth)} else { text(x = 59.5, y = min.sz + startmin2 + sh * sh.fac * seq(5, 0, by = -1), labels = formatC(1000 * MSEw, digits = 2, format = "f"), col = ifelse(MSEw == min(MSEw), "red", "black"), adj = 0.5, cex = texth)} anzahl2 <- apply(estw.temp, 2, function(x) sum(((x < yWerte.w[i, 1]) | (x > yWerte.w[i, 2])), na.rm = TRUE)) if(any(anzahl2)) { out.range.w[z2, ] <- c(i, l, anzahl2) z2 <- z2 + 1 } } # Eintragen der MSE- und Schätzer-Beschriftungen text(x = c(xlim2[1] - 0.3, 52.5), y = c(max.sz - startmax.sz, min.sz + startmin2) + sh * sh.fac * c(-0, 6), labels = paste("MSE * 1000:"), adj = 0) text(x = c(rep(xlim2[1] - 0.3, 6), rep(52.5, 6)), y = c(max.sz - startmax.sz - sh * sh.fac * seq(1, 6, by = 1), min.sz + startmin2 + sh * sh.fac * seq(5, 0, by = -1)), labels = rep(name, 2), cex = texth, adj = 0) ## Schreiben eines gemeinsamen Titels für beide Grafiken ############################################ if(Vplan[i, 1] != "Exp") { title(paste("Gruppe 0/Gruppe 1: ", Vplan[i,1], "(", Vplan[i,2], ", variierend)/", Vplan[i,1], "(", Vplan[i,2], ", ", round(as.numeric(Vplan[i,3]), 4), "), ", Vplan[i,4], "/", Vplan[i,5], " Beobachtungen", "\n Zensierungen: ", Vplan[i,6], ", Bindungen: ", Vplan[i,7], ", Störeffekt: ", Vplan[i,8], sep = ""), outer = TRUE) } else { title(paste("Gruppe 0/Gruppe 1: ", Vplan[i,1], "(variierend)/", Vplan[i,1], "(", round(as.numeric(Vplan[i,3]), 4), "), ", Vplan[i,4], "/", Vplan[i,5], " Beobachtungen", "\n Zensierungen: ", Vplan[i,6], ", Bindungen: ", Vplan[i,7], ", Störeffekt: ", Vplan[i,8], sep = ""), outer = TRUE) } dev.off() } save(out.range.lw, file = paste(file.erg, "/SimEstHR_Info_orBoxplotslw.RData", sep = "")) save(out.range.w, file = paste(file.erg, "/SimEstHR_Info_orBoxplotsw.RData", sep = "")) ### Zeichnung der beiden Grafiken für alle Szenarien in einer gemeinsamen Datei ################################### z1 <- 1 z2 <- 1 # Öffnen eines pdfs pdf(paste(file.erg, "/SimEstHR_Boxplots_Alle.pdf", sep = ""), width = 21/2.54, height = 29.7/2.54) for(i in 1:50) { ## Laden der Datei des entsprechenden Szenarios load(paste(file.erg, "/SimEstHR_Sz", i, ".RData", sep = "")) ## erste Grafik (für das logarithmierte Hazard-Ratio) ############################################### # Extrahierung der Schätzwerte aller sechs Schätzer est <- as.data.frame(lapply(sim.erg, "[", 1:6)) # allgemeine Einstellungen für das Layout der Grafiken par(mfrow = c(2, 1), lwd = 1.1, cex = 0.6, cex.axis = 1.4, cex.lab = 1.4, cex.main = 1.4, mar = c(4.5, 5, 0.5, 1.5), oma = c(0.1, 0.1, 3.5, 0.1)) min.sz <- min.1[which(Samp.Size[, 1] == Vplan$Stipro.G0[i])] max.sz <- max.1[which(Samp.Size[, 1] == Vplan$Stipro.G0[i])] startmax.sz <- startmax1[which(Samp.Size[, 1] == Vplan$Stipro.G0[i])] # Zeichnen der parallelen Boxplots boxplot(est, at = unlist(lapply(log.haz.rat, function(x) x + fuer.at * 0.25)), xlim = xlim1, ylim = c(min.sz, max.sz), las = 1, pars = list(boxwex = 0.017, staplewex = 0.5), col = mycols, xaxt = "n", cex = 0.2, xlab = expression(Wahres ~~ ln(omega)), ylab = expression(Geschätztes ~~ ln(omega))) axis(1, at = log.haz.rat, labels = round(log.haz.rat, 2)) for(j in seq(along = haz.rat)) { lines(c(log.haz.rat[j] + fuer.at[1] * 0.25 - 0.02, log.haz.rat[j] + fuer.at[6] * 0.25 + 0.02), rep(log.haz.rat[j], 2)) } # Schreiben einer Legende für die Zuordnung der Schätzer zu den einzelnen Boxplots if(Vplan$Bdg[i] == "keine") { legend("topleft", horiz = TRUE, legend = name, col = mycols, lwd = 1.8, lty = "solid", cex = 1.2) } else { legend("topleft", horiz = TRUE, legend = name.B, col = mycols, lwd = 1.8, lty = "solid", cex = 1.2) } # Eintragen der gruppenweisen Zensierungsanteile in die Grafik, wenn Zensierungen vorliegen texth <- 0.95 sh <- strheight(expression(S[W]), cex = texth) sh.fac <- 1.5 if(Vplan$Zens[i] != "keine") { n0 <- as.numeric(Vplan$Stipro.G0[i]); n1 <- as.numeric(Vplan$Stipro.G1[i]) zens0 <- (n0 - unlist(lapply(lapply(sim.erg, "[", 9), colMeans))) / n0 * 100 zens1 <- (n1 - unlist(lapply(lapply(sim.erg, "[", 10), colMeans))) / n1 * 100 text(x = c(xlim1[1] - 0.05, log.haz.rat[5]), y = c(max.sz - startmax.sz, min.sz + startmin1) + sh * sh.fac * c(-0, 8), labels = "Zens(%): ", cex = texth, adj = 0) text(x = log.haz.rat, y = rep(c(max.sz - startmax.sz, min.sz + startmin1) + sh * sh.fac * c(-0, 8), c(5, 4)), labels = paste(formatC(zens0, digits = 0, format = "f"), "/", formatC(zens1, digits = 0, format = "f"), sep = ""), cex = texth) } # Eintragen des gruppenübergreifenden Bindungssanteils in die Grafik, wenn Bindungen vorliegen if(Vplan$Bdg[i] != "keine") { n <- as.numeric(Vplan$Stipro.G0[i]) + as.numeric(Vplan$Stipro.G1[i]) bd01 <- numeric(length(haz.rat)) for(l in seq(along = haz.rat)) { sim.erg.l <- sim.erg[[l]] bd01[l] <- mean(sim.erg.l[, 13] / sim.erg.l[, 8]) * 100 } text(x = c(xlim1[1] - 0.05, log.haz.rat[5]), y = c(max.sz - startmax.sz, min.sz + startmin1) + sh * sh.fac * c(-1, 7), labels = "Bdg(%): ", cex = texth, adj = 0) text(x = log.haz.rat, y = rep(c(max.sz - startmax.sz, min.sz + startmin1) + sh * sh.fac * c(-1, 7), c(5, 4)), labels = paste(formatC(bd01, digits = 0, format = "f"), sep = ""), cex = texth) } # Eintragen der MSEs in die Grafik für die einzelnen logarithmierten Hazard-Ratios und Schätzer sowie # Speicherung der Anzahl d. Punkte außerhalb des Plotbereichs pro Schätzer für die einzelnen log. Hazard-Ratios yWerte.lw[i, ] <- par("usr")[c(3, 4)] for(l in seq(along = haz.rat)) { est.temp <- sim.erg[[l]][,1:6] MSE <- apply(est.temp, 2, function(x) var(x, na.rm = TRUE)) + (apply(est.temp, 2, function(x) mean(x, na.rm = TRUE)) - log.haz.rat[l])^2 max.MSE <- which(MSE == max(MSE)) if(l <= 5) { text(x = log.haz.rat[l], y = max.sz - startmax.sz - sh * sh.fac * seq(3, 8, by = 1), labels = formatC(1000 * MSE, digits = 2, format = "f"), col = ifelse(MSE == min(MSE), "red", "black"), adj = 0.5, cex = texth)} else { text(x = log.haz.rat[l], y = min.sz + startmin1 + sh * sh.fac * seq(5, 0, by = -1), labels = formatC(1000 * MSE, digits = 2, format = "f"), col = ifelse(MSE == min(MSE), "red", "black"), adj = 0.5, cex = texth)} anzahl1 <- apply(est.temp, 2, function(x) sum(((x < yWerte.lw[i, 1]) | (x > yWerte.lw[i, 2])), na.rm = TRUE)) if(any(anzahl1)) { out.range.lw[z1, ] <- c(i, l, anzahl1) z1 <- z1 + 1 } } # Eintragen der MSE- und Schätzer-Beschriftungen text(x = c(xlim1[1] - 0.05, log.haz.rat[5]), y = c(max.sz - startmax.sz, min.sz + startmin1) + sh * sh.fac * c(-2, 6), labels = paste("MSE * 1000:"), adj = 0) text(x = c(rep(xlim1[1] - 0.05, 6), rep(log.haz.rat[5], 6)), y = c(max.sz - startmax.sz - sh * sh.fac * seq(3, 8, by = 1), min.sz + startmin1 + sh * sh.fac * seq(5, 0, by = -1)), labels = rep(name, 2), cex = texth, adj = 0) ## zweite Grafik (für das nicht-logarithmierte Hazard-Ratio) ######################################## # Extrahierung Schätzwerte aller sechs Schätzer estw <- exp(est) # allgemeine Einstellungen für das Layout der Grafiken par(mar = c(4.5, 5, 2.5, 1.5)) min.sz <- min.2[which(Samp.Size[, 1] == Vplan$Stipro.G0[i])] max.sz <- max.2[which(Samp.Size[, 1] == Vplan$Stipro.G0[i])] startmax.sz <- startmax2[which(Samp.Size[, 1] == Vplan$Stipro.G0[i])] # Zeichnen der parallelen Boxplots boxplot(estw, at = xaxis.2, xlim = c(0.9, 62.1), ylim = c(min.sz, max.sz), las = 1, pars = list(boxwex = 0.7, staplewex = 0.5), col = mycols, xaxt = "n", cex = 0.2, xlab = expression(Wahres ~~ omega), ylab = expression(Geschätztes ~~ omega)) axis(1, at = seq(3.5, 59.5, by = 7), labels = round(haz.rat, 2)) for(k in seq(along = haz.rat)) { lines(c(seq(1, 57, by = 7)[k] - 0.3, seq(6, 62, by = 7)[k]) + 0.3, rep(haz.rat[k], 2)) } # Schreiben einer Legende für die Zuordnung der Schätzer zu den einzelnen Boxplots if(Vplan$Bdg[i] == "keine") { legend("topleft", horiz = TRUE, legend = name.exp, col = mycols, lwd = 1.8, lty = "solid", cex = 1.2) } else { legend("topleft", horiz = TRUE, legend = name.expB, col = mycols, lwd = 1.8, lty = "solid", cex = 1.2) } # Eintragen der MSEs in die Grafik für die einzelnen Hazard-Ratios und Schätzer sowie # Speicherung der Anzahl d. Punkte außerhalb des Plotbereichs pro Schätzer für die einzelnen Hazard-Ratios sh <- strheight(expression(S[W]), cex = texth) yWerte.w[i, ] <- par("usr")[c(3, 4)] for(l in seq(along = haz.rat)) { estw.temp <- exp(sim.erg[[l]][, 1:6]) MSEw <- apply(estw.temp, 2, function(x) var(x, na.rm = TRUE)) + (apply(estw.temp, 2, function(x) mean(x, na.rm = TRUE)) - haz.rat[l])^2 if(l <= 8) { text(x = xaxis.2.re[l], y = max.sz - startmax.sz - sh * sh.fac * seq(1, 6, by = 1), labels = formatC(1000 * MSEw, digits = 2, format = "f"), col = ifelse(MSEw == min(MSEw), "red", "black"), adj = 0.5, cex = texth)} else { text(x = 59.5, y = min.sz + startmin2 + sh * sh.fac * seq(5, 0, by = -1), labels = formatC(1000 * MSEw, digits = 2, format = "f"), col = ifelse(MSEw == min(MSEw), "red", "black"), adj = 0.5, cex = texth)} anzahl2 <- apply(estw.temp, 2, function(x) sum(((x < yWerte.w[i, 1]) | (x > yWerte.w[i, 2])), na.rm = TRUE)) if(any(anzahl2)) { out.range.w[z2, ] <- c(i, l, anzahl2) z2 <- z2 + 1 } } # Eintragen der MSE- und Schätzer-Beschriftungen text(x = c(xlim2[1] - 0.3, 52.5), y = c(max.sz - startmax.sz, min.sz + startmin2) + sh * sh.fac * c(-0, 6), labels = paste("MSE * 1000:"), adj = 0) text(x = c(rep(xlim2[1] - 0.3, 6), rep(52.5, 6)), y = c(max.sz - startmax.sz - sh * sh.fac * seq(1, 6, by = 1), min.sz + startmin2 + sh * sh.fac * seq(5, 0, by = -1)), labels = rep(name, 2), cex = texth, adj = 0) ## Schreiben eines gemeinsamen Titels für beide Grafiken ############################################ if(Vplan[i, 1] != "Exp") { title(paste("Gruppe 0/Gruppe 1: ", Vplan[i,1], "(", Vplan[i,2], ", variierend)/", Vplan[i,1], "(", Vplan[i,2], ", ", round(as.numeric(Vplan[i,3]), 4), "), ", Vplan[i,4], "/", Vplan[i,5], " Beobachtungen", "\n Zensierungen: ", Vplan[i,6], ", Bindungen: ", Vplan[i,7], ", Störeffekt: ", Vplan[i,8], sep = ""), outer = TRUE) } else { title(paste("Gruppe 0/Gruppe 1: ", Vplan[i,1], "(variierend)/", Vplan[i,1], "(", round(as.numeric(Vplan[i,3]), 4), "), ", Vplan[i,4], "/", Vplan[i,5], " Beobachtungen", "\n Zensierungen: ", Vplan[i,6], ", Bindungen: ", Vplan[i,7], ", Störeffekt: ", Vplan[i,8], sep = ""), outer = TRUE) } } dev.off() }