Nothing
R/report.r
In eatRep: Educational Assessment Tools for Replication Methods
Defines functions
reduceDoubleN
compare_point_estimates
seCorrect.rep_se_correction
recursive_filter
separate_jk2
seCorrect.wec_se_correction
seCorrect.old
seCorrect
report <- function ( repFunOut, trendDiffs = FALSE, add=list(), exclude = c("NcasesValid", "var", "sampleSize"), printGlm = FALSE,
round = TRUE, digits = 3, printDeviance = FALSE) {
if(is.null(repFunOut)) {return(NULL)}
### vorab: alte 'dG'-Funktion zum Anzeigen der Regressionsergebnisse implementieren
if ( length(grep("glm", as.character(repFunOut[["resT"]][[1]][1,"modus"]))) ==1 ) {
if ( printGlm == TRUE ) { dG(repFunOut, digits = digits, printDeviance = printDeviance, add = add ) }
}
### 1. Input extrahieren: diese Variablen dann spaeter an Einzelfunktionen weitergeben!
jk2 <- repFunOut[["resT"]]
tv <- repFunOut[["allNam"]][["trend"]]
cols <- c("group", "depVar", "modus", "parameter")
grpv <- setdiff(setdiff(colnames(jk2[[1]]), cols), c("comparison", "coefficient", "value", tv))
grp_by <- repFunOut[["allNam"]][["group.differences.by"]]
cl_diffs <- repFunOut[["allNam"]][["cross.differences"]]
funs <- c("mean", "table", "quantile", "glm", "lmer")
fun <- funs [ which( unlist(lapply(funs, FUN = function ( f ) { length(grep(f, jk2[[1]][1,"modus"]))})) > 0) ]
### 2. cross-level diffs bestimmen: ueberschreibt bzw. erweitert das Objekt 'jk2' ... Achtung: sind nur fuer "mean" oder "table" erlaubt
if ( is.list(cl_diffs) ) {
jk2 <- lapply(jk2, FUN = function (df) {computeCrossLevel (df, cols=cols, grpv = grpv, fun = fun, cl_diffs = cl_diffs, comp_type = "crossDiff")})
jk2 <- lapply(jk2, FUN = function (df) { ### in spalte "comparison" 'crossDiff_of_groupDiff' eintragen, falls in Spalte "group" 3x ".vs." steht
spl <- strsplit(df[,"group"], ".vs.")
ind <- which(sapply(spl, length)==4)
if ( length(ind)>0) {
df[ind,"comparison"] <- "crossDiff_of_groupDiff" ### das mittlere ".vs." gross schreiben
df[ind,"group"] <- unlist(lapply(spl[ind], FUN = function ( z ) {paste(z[1], ".vs.", z[2], ".VS.", z[3], ".vs.", z[4], sep="")}))
}
return(df)})
}
### 2.b) SE correction durchfuehren (siehe Paper Weirich & Hecht)
if(!is.null(repFunOut[["SE_correction"]]) && !is.null(repFunOut[["SE_correction"]][[1]])) {
### checks, ob Vergleiche dabei, fuer die keine Korrektur verfuehgbar ist
if(length(which(jk2[[1]][["comparison"]] == "crossDiff_of_groupDiff")) > 0 ) {
warning("Standard error correction for 'crossDiff_of_groupDiff' is currently not supported.")
}
mult_hierarchy <- any(unlist(lapply(repFunOut$allNam$cross.differences, function(x) x[2] - x[1] != 1)))
if(mult_hierarchy) warning("Standard error correction for crossDifferences across multiple hierarchy levels is currently not supported.")
### correction durchfuehren
jk2 <- lapply(jk2, function(jk2_single) { seCorrect(SE_correction = repFunOut[["SE_correction"]], jk2 = jk2_single, grpv = grpv) })
}
### 3. Trend bestimmen
if ( !is.null(tv) ) {
jk2 <- computeTrend(jk2 = jk2, repFunOut = repFunOut, tv = tv, fun = fun)
} else {
jk2 <- jk2[[1]]
}
### 4. Trend-Differences bestimmen ... Achtung: sind nur fuer "mean" oder "table" erlaubt
if ( !is.null(tv) && trendDiffs ) {
jk2 <- computeTrendDiffs(jk2 = jk2, grpv = grpv, tv = tv, grp_by = grp_by, fun = fun, cl_diffs = cl_diffs)
}
### 5. 'add' ergaenzen, falls gewuenscht
if ( length(add)>0) {
if(!all(nchar(names(add))>0)) { stop("'add' must be named.")} ### necessary checks
if(length(names(add)) != length(unique(names(add)))) { stop("Duplicated names of 'add' are not allowed.")}
if(!all(sapply(add, length) == 1)) {stop("All elements of 'add' must be of length 1.")}
if(!all(sapply(add, class) == "character")) {stop("All elements of 'add' must be of class 'character'.")}
dopp<- names(add) %in% colnames(jk2)
ind <- which(dopp==TRUE)
if ( length( ind ) > 0 ) {stop(paste0("Following names of 'add' are not allowed: '",paste(names(add)[ind], collapse = "', '"), "'."))}
for ( u in names(add)) {jk2[,u] <- add[[u]]}
}
### 6. reshapen
spltVar <- c("coefficient", tv) ### split-Variable
if ( length(exclude)>0) {
weg <- which(jk2[,"parameter"] %in% exclude)
if ( length(weg)>0) {
jk2 <- jk2[-weg,]
}
} ### was muss in die Spalten? das haengt davon ab, ob es einen Trend gibt
frml <- as.formula(paste0("... ~ ", paste(spltVar,collapse=" + ") ) )
### Hotfix: wenn repTable ueber wiederholten Aufrufen von repMean gewrappt wurde, stehen die Ns mehrmals drin, naemlich fuer
### jede Indikatorvariable einer mehrstufigen Faktorvariable separat. Sie sind aber immer gleich. Durch das Mehrmalsdrinstehen
### misslingt das reshapen, deshalb muessen sie jetzt raus
if ( fun == "table") {jk2 <- reduceDoubleN(jk2)}
### Hotfix: damit fuer glm im Output die Koeffizienten immer zuerst und R2, R2nagel, Nvalid imer dahinter stehen, werden die jetzt voruebergehen umbenannt (hinterher wieder zurueck benannt)
if ( fun == "glm") {jk2[,"parameter"] <- car::recode(jk2[,"parameter"], "'Nvalid'='zzzzNvalid'; 'R2'='zzzzR2'; 'R2nagel'='zzzzR2nagel'") }
jk2wide <- reshape2::dcast(data = jk2, formula = frml, value.var = "value")
if ( fun == "glm") {jk2wide[,"parameter"] <- car::recode(jk2wide[,"parameter"], "'zzzzNvalid'='Nvalid'; 'zzzzR2'='R2'; 'zzzzR2nagel'='R2nagel'") }
### runden, falls gewuenscht
if ( isTRUE(round)) {
jk2wide <- eatTools::roundDF(jk2wide, digits = digits)
}
return(jk2wide)}
addSig <- function ( dat , groupCols = NULL , allNam = NULL ) {
if(is.null(groupCols)) {groupCols <- c("group", "parameter")}
dat <- do.call("rbind", by ( data = dat, INDICES = dat[,groupCols], FUN = function ( x ) {
z <- x[which(x[,"coefficient"] %in% c("est", "se")),]
if ( nrow(z) > 2) {cat("Fehler. x muss maximal 2 zeilen haben.\n")}
if ( nrow(z) == 2 ) {
y <- z[1,] ### dazu relevante spalten identifizieren, nach denen gesplittet werden muss
y[["coefficient"]] <- "p" ### erste Zeile von x duplizieren und relevante Werte ersetzen
y[["value"]] <- 2*pnorm(abs(z[which(z[,"coefficient"] == "est"),"value"] / z[which(z[,"coefficient"] == "se"),"value"]), lower.tail=FALSE)
x <- rbind ( x, y) ### Achtung: Signifikanzwert wird hier noch nach numerisch transformiert, muss zurueckgewandelt werden
}
return(x)})) ### untere Zeile: wenn 'table' ueber 'jk2.mean' gewrappt wurde, muessen hier die parameterbezeichnungen geaendert werden
return(dat)}
seCorrect <- function( SE_correction, jk2, grpv ) {
# stop("SE correction has not been implemented yet. Use crossDiffSE = 'old'.")
UseMethod("seCorrect")
}
seCorrect.old <- function( SE_correction, jk2, grpv ) {
jk2
}
## an der falschen Stelle! muss vor Trends passieren!
seCorrect.wec_se_correction <- function( SE_correction, jk2, grpv ) {
sep_jk2 <- separate_jk2(jk2 = jk2)
cross_diff <- sep_jk2[["cross_diff"]]
for(i in seq_along(SE_correction)) {
# debug 3 grp-vars: i <- 4
output <- SE_correction[[i]][["resT"]][[sep_jk2[["year"]]]]
rows <- length(SE_correction[[i]][["vgl"]][["groups.divided.by"]])
single_grpv <- SE_correction[[i]][["focGrp"]]
# Nicht output reporten, da spalten sonst je nach Trend unterschiedlich heissen
SEs <- output[!output$parameter %in% c("(Intercept)", "Nvalid", "R2"), c("parameter", "value", "coefficient")]
SEs[, "parameter"] <- gsub(paste0("^", single_grpv), "", SEs[, "parameter"])
SEs <- as.data.frame(tidyr::pivot_wider(SEs, names_from = "coefficient", values_from = "value"))
#if(is.data.frame(SE_correction[[i]]$refGrp) && identical(SE_correction[[i]]$refGrp$groupValue, "LandA")) browser()
for(param in SEs[["parameter"]]) {
esc_param <- Hmisc::escapeRegex(param)
if(identical(SE_correction[[i]]$refGrp, "all")) { ## if reference level is the whole group
grp_regexp <- paste0("^", esc_param, "\\.vs")
compare_point_estimates(old_est = cross_diff[cross_diff$parameter == "mean" & grepl(grp_regexp, cross_diff$group) & cross_diff$coefficient == "est", "value"],
new_est = SEs[SEs[, "parameter"] == param, "est"],
param = param)
# funktioniert nur, da keine correction ueber mehrere Hierarchieebenen hinweg funktioniert!
cross_diff[cross_diff$parameter == "mean" & grepl(grp_regexp, cross_diff$group) & cross_diff$coefficient == "se",
"value"] <- SEs[SEs[, "parameter"] == param, "se"]
cross_diff[cross_diff$parameter == "mean" & grepl(grp_regexp, cross_diff$group) & cross_diff$coefficient == "p",
"value"] <- SEs[SEs[, "parameter"] == param, "p"]
} else { ### if reference level is a subgroup
# create variable to select right rows in jk2 output
#param_finder <- cross_diff$group
#param_finder <- sapply(strsplit(param_finder, ".vs."), function(x) x[1])
## complicated to find param match, because the factor string of another level can contain the string of the current level!
param_selector <- paste0("^", esc_param, "\\.|", "^", esc_param, "_|", "_", esc_param, "\\.|", "_", esc_param, "_")
col_names <- SE_correction[[i]]$refGrp[, "groupName"]
col_levels <- SE_correction[[i]]$refGrp[, "groupValue"]
#if(identical(col_levels, c("female", "TRUE"))) browser()
# filter relevant rows for flexibel number of filter variables
filt_var <- recursive_filter(df = cross_diff, vars = col_names, var_levels = col_levels)
compare_point_estimates(old_est = cross_diff[which(filt_var & grepl(param_selector, cross_diff$group) & cross_diff$coefficient == "est"), "value"],
new_est = SEs[SEs[, "parameter"] == param, "est"],
param = param)
cross_diff[which(filt_var & grepl(param_selector, cross_diff$group) & cross_diff$coefficient == "se"),
"value"] <- SEs[SEs[, "parameter"] == param, "se"]
cross_diff[which(filt_var & grepl(param_selector, cross_diff$group) & cross_diff$coefficient == "p"),
"value"] <- SEs[SEs[, "parameter"] == param, "p"]
}
}
}
# change modus in reporting output for WEC
cross_diff$modus <- unique(SE_correction[[1]][[1]][[1]]$modus)
## mehr Fragen
# _ in Faktorlevels erlaubt?
# next: trend, groupdiffs checken, ob hier SEs korrekt beruecksichtigt werden
# was passiert eig. wenn faktorlevel verschiedener Variablen gleich heissen?
rbind(sep_jk2[["no_cross_diff"]], cross_diff)
}
## Prepare jk2 object for se_correction
separate_jk2 <- function(jk2) {
## Trend or no trend
year <- as.character(unique(jk2[["year"]]))
if(length(year) == 0) year <- 1
## dissect results so far
no_cross_diff <- jk2[is.na(jk2$comparison) | jk2$comparison != "crossDiff" | jk2$parameter != "mean", ]
cross_diff <- jk2[which(jk2$comparison == "crossDiff" & jk2$parameter == "mean"), ]
stopifnot(identical(nrow(no_cross_diff) + nrow(cross_diff), nrow(jk2)))
stopifnot(nrow(cross_diff) > 0)
list(cross_diff = cross_diff, no_cross_diff = no_cross_diff, year = year)
}
recursive_filter <- function(df, vars, var_levels) {
stopifnot(length(vars) == length(var_levels))
filt_var <- rep(TRUE, nrow(df))
for(i in seq_along(vars)) {
filt_var <- filt_var & df[[vars[i]]] == var_levels[[i]]
}
# output: vector with T/F
filt_var
}
seCorrect.rep_se_correction <- function( SE_correction, jk2, grpv ) {
sep_jk2 <- separate_jk2(jk2 = jk2)
cross_diff <- sep_jk2[["cross_diff"]]
no_cross_diff <- sep_jk2[["no_cross_diff"]]
for(i in seq_along(SE_correction)) {
output <- SE_correction[[i]][["resT"]][[sep_jk2[["year"]]]]
rows <- length(SE_correction[[i]][["vgl"]][["groups.divided.by"]])
single_grpv <- SE_correction[[i]][["focGrp"]]
for(param in output[["group"]]) {
compare_point_estimates(old_est = cross_diff[cross_diff$group == param & cross_diff$coefficient == "est", "value"],
new_est = output[output$group == param & output$coefficient == "est", "value"],
param = param)
cross_diff[cross_diff$group == param & cross_diff$coefficient == "se", "value"] <- output[output$group == param & output$coefficient == "se", "value"]
cross_diff[cross_diff$group == param & cross_diff$coefficient == "p", "value"] <- output[output$group == param & output$coefficient == "p", "value"]
}
}
rbind(no_cross_diff, cross_diff)
}
compare_point_estimates <- function(old_est, new_est, param) {
if(abs(abs(old_est) - abs(new_est)) >= 1e-10) {
#browser()
warning("Difference in point estimate of cross level difference for comparison ", param, ": ", round(abs(old_est) - abs(new_est), digits = 3))
}
return()
}
### Hier kann man nicht subset() nehmen, weil es sonst Warnungen beim Paketebauen gibt,
### da subset() die Spaltennamen unquoted haben will
reduceDoubleN <- function(jk2){
cases <- unique(jk2[which(jk2[,"parameter"] == "Ncases"),]) ### cases <- unique(subset(jk2, parameter=="Ncases"))
jk2 <- rbind(jk2[which(jk2[,"parameter"] != "Ncases"),], cases) ### jk2 <- rbind(subset(jk2,parameter != "Ncases"), cases)
return(jk2)}
#reshapeReport <- function(inp, tv, fun, repFunOut, target, wholeGroupName) {
### testweise initialisieren
# inpR <- inp ### untere zeile: levels der Trendvariablen
# lev <- paste(as.character(unique(do.call("rbind", repFunOut[["resT"]])[,repFunOut[["allNam"]][["trend"]]])), collapse="|")
### reshapen nur machen, wenn Parameterspalte mehr als einen Eintrag hat ... das ist bspw. nicht der Fall, wenn man repTable() mit dichotomer AV (Regelstandard erreicht) macht
# if (length(unique(inp[,"parameter"])) > 1) {
### fuer trend
# if (!is.null(tv) ) {
### fuer mittelwerte mit trend
# if (fun == "mean") {
# inpR <- data.frame(tidyr::pivot_wider(inp, names_from = "parameter", values_from = grep(lev, colnames(inp), value=TRUE)), stringsAsFactors=FALSE)
# }
### fuer haeufigkeitsverteilungen mit trend
# if (fun == "table") {
# inpR <- data.frame(tidyr::pivot_wider(inp, names_from = "parameter", values_from = grep(lev, colnames(inp), value=TRUE)), stringsAsFactors=FALSE)
# }
# } else {
### fuer mittelwerte ohne trend
# if (fun == "mean") {
# inpR <- data.frame(tidyr::pivot_wider(inp, names_from = "parameter", values_from = intersect(c("es", "est", "p", "se"), colnames(inp))), stringsAsFactors=FALSE)
# }
### fuer haeufigkeitsverteilungen ohne trend
# if (fun == "table") {
# inpR <- data.frame(tidyr::pivot_wider(inp, names_from = "parameter", values_from = intersect(c("est", "p", "se"), colnames(inp))), stringsAsFactors=FALSE)
# }
# }
# }
### fuer allen Output: spalten, die nur NA haben, loeschen
# inpR <- inpR[, colSums(is.na(inpR)) != nrow(inpR)]
### crossDiff-Benennung anpassen (Mail Juan Jose, 27.04.2022, 9.37 Uhr). Das passiert auch dann, wenn nicht gereshapet wird
# if ("comparison" %in% colnames(inpR)) {
# if( any(c("crossDiff", "trendDiff_cross") %in% inpR[,"comparison"])) {
# rows <- intersect(intersect(grep("wholeGroup$", inpR[,"group"]), setdiff(1:nrow(inpR), grep("_",inpR[,"group"]))), which(inpR[,"comparison"] %in% c("crossDiff", "trendDiff_cross")))
# if (length(rows)>0) {
# inpR[rows,"group"] <- eatTools::removePattern(inpR[rows,"group"], "\\.vs\\.wholeGroup$")
# }
# }
# }
### Alternativ: Mail Kristoph, 12.05.2022, 18.40 Uhr (geht wahrscheinlich nur fuer eine Gruppierungsvariable)
# if (fun == "table" && length(unique(inpR[,"group"])) < nrow(inpR) && target == "BT2021.table" && length(repFunOut[["allNam"]][["group"]])>0 ) {
# prz <- grep("^est|^se", colnames(inpR), value=TRUE)
# if(length(prz)<1) {warning("Cannot found colums with percentage values. Skip value transformation.")} else {
# for ( i in prz) { inpR[,i] <- 100*inpR[,i]}
# }
### Spalten identifizieren, die fuers reshapen entfernt werden muessen, und zwar sowohl fuer mit als auch fuer ohne trend
# weg <- grep(repFunOut[["allNam"]][["group"]], colnames(inpR))
# stopifnot(length(weg)>0)
### erstmal ohne Trend
# if (is.null(tv) ) {
# inpR <- data.frame(tidyr::pivot_wider(inpR[,-weg], names_from = "comparison", values_from = intersect(c("es", "est", "p", "se"), colnames(inpR))), stringsAsFactors=FALSE)
# } else { ### mit Trend
# inpR <- data.frame(tidyr::pivot_wider(inpR[,-weg], names_from = "comparison", values_from = grep(lev, colnames(inp), value=TRUE)), stringsAsFactors=FALSE)
# }
# inpR[,"group"] <- car::recode(inpR[,"group"], paste0("'wholeGroup'='",wholeGroupName,"'"))
# colnames(inpR) <- eatTools::removePattern(colnames(inpR), pattern = "_NA$")
# }
# inpR <- inpR[, colSums(is.na(inpR)) != nrow(inpR)] ### nochmal spalten loeschen mit nur NA
# return(inpR)}
Try the eatRep package in your browser
Any scripts or data that you put into this service are public.
eatRep documentation built on March 31, 2023, 7:15 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions reduceDoubleN compare_point_estimates seCorrect.rep_se_correction recursive_filter separate_jk2 seCorrect.wec_se_correction seCorrect.old seCorrect
report <- function ( repFunOut, trendDiffs = FALSE, add=list(), exclude = c("NcasesValid", "var", "sampleSize"), printGlm = FALSE,
round = TRUE, digits = 3, printDeviance = FALSE) {
if(is.null(repFunOut)) {return(NULL)}
### vorab: alte 'dG'-Funktion zum Anzeigen der Regressionsergebnisse implementieren
if ( length(grep("glm", as.character(repFunOut[["resT"]][[1]][1,"modus"]))) ==1 ) {
if ( printGlm == TRUE ) { dG(repFunOut, digits = digits, printDeviance = printDeviance, add = add ) }
}
### 1. Input extrahieren: diese Variablen dann spaeter an Einzelfunktionen weitergeben!
jk2 <- repFunOut[["resT"]]
tv <- repFunOut[["allNam"]][["trend"]]
cols <- c("group", "depVar", "modus", "parameter")
grpv <- setdiff(setdiff(colnames(jk2[[1]]), cols), c("comparison", "coefficient", "value", tv))
grp_by <- repFunOut[["allNam"]][["group.differences.by"]]
cl_diffs <- repFunOut[["allNam"]][["cross.differences"]]
funs <- c("mean", "table", "quantile", "glm", "lmer")
fun <- funs [ which( unlist(lapply(funs, FUN = function ( f ) { length(grep(f, jk2[[1]][1,"modus"]))})) > 0) ]
### 2. cross-level diffs bestimmen: ueberschreibt bzw. erweitert das Objekt 'jk2' ... Achtung: sind nur fuer "mean" oder "table" erlaubt
if ( is.list(cl_diffs) ) {
jk2 <- lapply(jk2, FUN = function (df) {computeCrossLevel (df, cols=cols, grpv = grpv, fun = fun, cl_diffs = cl_diffs, comp_type = "crossDiff")})
jk2 <- lapply(jk2, FUN = function (df) { ### in spalte "comparison" 'crossDiff_of_groupDiff' eintragen, falls in Spalte "group" 3x ".vs." steht
spl <- strsplit(df[,"group"], ".vs.")
ind <- which(sapply(spl, length)==4)
if ( length(ind)>0) {
df[ind,"comparison"] <- "crossDiff_of_groupDiff" ### das mittlere ".vs." gross schreiben
df[ind,"group"] <- unlist(lapply(spl[ind], FUN = function ( z ) {paste(z[1], ".vs.", z[2], ".VS.", z[3], ".vs.", z[4], sep="")}))
}
return(df)})
}
### 2.b) SE correction durchfuehren (siehe Paper Weirich & Hecht)
if(!is.null(repFunOut[["SE_correction"]]) && !is.null(repFunOut[["SE_correction"]][[1]])) {
### checks, ob Vergleiche dabei, fuer die keine Korrektur verfuehgbar ist
if(length(which(jk2[[1]][["comparison"]] == "crossDiff_of_groupDiff")) > 0 ) {
warning("Standard error correction for 'crossDiff_of_groupDiff' is currently not supported.")
}
mult_hierarchy <- any(unlist(lapply(repFunOut$allNam$cross.differences, function(x) x[2] - x[1] != 1)))
if(mult_hierarchy) warning("Standard error correction for crossDifferences across multiple hierarchy levels is currently not supported.")
### correction durchfuehren
jk2 <- lapply(jk2, function(jk2_single) { seCorrect(SE_correction = repFunOut[["SE_correction"]], jk2 = jk2_single, grpv = grpv) })
}
### 3. Trend bestimmen
if ( !is.null(tv) ) {
jk2 <- computeTrend(jk2 = jk2, repFunOut = repFunOut, tv = tv, fun = fun)
} else {
jk2 <- jk2[[1]]
}
### 4. Trend-Differences bestimmen ... Achtung: sind nur fuer "mean" oder "table" erlaubt
if ( !is.null(tv) && trendDiffs ) {
jk2 <- computeTrendDiffs(jk2 = jk2, grpv = grpv, tv = tv, grp_by = grp_by, fun = fun, cl_diffs = cl_diffs)
}
### 5. 'add' ergaenzen, falls gewuenscht
if ( length(add)>0) {
if(!all(nchar(names(add))>0)) { stop("'add' must be named.")} ### necessary checks
if(length(names(add)) != length(unique(names(add)))) { stop("Duplicated names of 'add' are not allowed.")}
if(!all(sapply(add, length) == 1)) {stop("All elements of 'add' must be of length 1.")}
if(!all(sapply(add, class) == "character")) {stop("All elements of 'add' must be of class 'character'.")}
dopp<- names(add) %in% colnames(jk2)
ind <- which(dopp==TRUE)
if ( length( ind ) > 0 ) {stop(paste0("Following names of 'add' are not allowed: '",paste(names(add)[ind], collapse = "', '"), "'."))}
for ( u in names(add)) {jk2[,u] <- add[[u]]}
}
### 6. reshapen
spltVar <- c("coefficient", tv) ### split-Variable
if ( length(exclude)>0) {
weg <- which(jk2[,"parameter"] %in% exclude)
if ( length(weg)>0) {
jk2 <- jk2[-weg,]
}
} ### was muss in die Spalten? das haengt davon ab, ob es einen Trend gibt
frml <- as.formula(paste0("... ~ ", paste(spltVar,collapse=" + ") ) )
### Hotfix: wenn repTable ueber wiederholten Aufrufen von repMean gewrappt wurde, stehen die Ns mehrmals drin, naemlich fuer
### jede Indikatorvariable einer mehrstufigen Faktorvariable separat. Sie sind aber immer gleich. Durch das Mehrmalsdrinstehen
### misslingt das reshapen, deshalb muessen sie jetzt raus
if ( fun == "table") {jk2 <- reduceDoubleN(jk2)}
### Hotfix: damit fuer glm im Output die Koeffizienten immer zuerst und R2, R2nagel, Nvalid imer dahinter stehen, werden die jetzt voruebergehen umbenannt (hinterher wieder zurueck benannt)
if ( fun == "glm") {jk2[,"parameter"] <- car::recode(jk2[,"parameter"], "'Nvalid'='zzzzNvalid'; 'R2'='zzzzR2'; 'R2nagel'='zzzzR2nagel'") }
jk2wide <- reshape2::dcast(data = jk2, formula = frml, value.var = "value")
if ( fun == "glm") {jk2wide[,"parameter"] <- car::recode(jk2wide[,"parameter"], "'zzzzNvalid'='Nvalid'; 'zzzzR2'='R2'; 'zzzzR2nagel'='R2nagel'") }
### runden, falls gewuenscht
if ( isTRUE(round)) {
jk2wide <- eatTools::roundDF(jk2wide, digits = digits)
}
return(jk2wide)}
addSig <- function ( dat , groupCols = NULL , allNam = NULL ) {
if(is.null(groupCols)) {groupCols <- c("group", "parameter")}
dat <- do.call("rbind", by ( data = dat, INDICES = dat[,groupCols], FUN = function ( x ) {
z <- x[which(x[,"coefficient"] %in% c("est", "se")),]
if ( nrow(z) > 2) {cat("Fehler. x muss maximal 2 zeilen haben.\n")}
if ( nrow(z) == 2 ) {
y <- z[1,] ### dazu relevante spalten identifizieren, nach denen gesplittet werden muss
y[["coefficient"]] <- "p" ### erste Zeile von x duplizieren und relevante Werte ersetzen
y[["value"]] <- 2*pnorm(abs(z[which(z[,"coefficient"] == "est"),"value"] / z[which(z[,"coefficient"] == "se"),"value"]), lower.tail=FALSE)
x <- rbind ( x, y) ### Achtung: Signifikanzwert wird hier noch nach numerisch transformiert, muss zurueckgewandelt werden
}
return(x)})) ### untere Zeile: wenn 'table' ueber 'jk2.mean' gewrappt wurde, muessen hier die parameterbezeichnungen geaendert werden
return(dat)}
seCorrect <- function( SE_correction, jk2, grpv ) {
# stop("SE correction has not been implemented yet. Use crossDiffSE = 'old'.")
UseMethod("seCorrect")
}
seCorrect.old <- function( SE_correction, jk2, grpv ) {
jk2
}
## an der falschen Stelle! muss vor Trends passieren!
seCorrect.wec_se_correction <- function( SE_correction, jk2, grpv ) {
sep_jk2 <- separate_jk2(jk2 = jk2)
cross_diff <- sep_jk2[["cross_diff"]]
for(i in seq_along(SE_correction)) {
# debug 3 grp-vars: i <- 4
output <- SE_correction[[i]][["resT"]][[sep_jk2[["year"]]]]
rows <- length(SE_correction[[i]][["vgl"]][["groups.divided.by"]])
single_grpv <- SE_correction[[i]][["focGrp"]]
# Nicht output reporten, da spalten sonst je nach Trend unterschiedlich heissen
SEs <- output[!output$parameter %in% c("(Intercept)", "Nvalid", "R2"), c("parameter", "value", "coefficient")]
SEs[, "parameter"] <- gsub(paste0("^", single_grpv), "", SEs[, "parameter"])
SEs <- as.data.frame(tidyr::pivot_wider(SEs, names_from = "coefficient", values_from = "value"))
#if(is.data.frame(SE_correction[[i]]$refGrp) && identical(SE_correction[[i]]$refGrp$groupValue, "LandA")) browser()
for(param in SEs[["parameter"]]) {
esc_param <- Hmisc::escapeRegex(param)
if(identical(SE_correction[[i]]$refGrp, "all")) { ## if reference level is the whole group
grp_regexp <- paste0("^", esc_param, "\\.vs")
compare_point_estimates(old_est = cross_diff[cross_diff$parameter == "mean" & grepl(grp_regexp, cross_diff$group) & cross_diff$coefficient == "est", "value"],
new_est = SEs[SEs[, "parameter"] == param, "est"],
param = param)
# funktioniert nur, da keine correction ueber mehrere Hierarchieebenen hinweg funktioniert!
cross_diff[cross_diff$parameter == "mean" & grepl(grp_regexp, cross_diff$group) & cross_diff$coefficient == "se",
"value"] <- SEs[SEs[, "parameter"] == param, "se"]
cross_diff[cross_diff$parameter == "mean" & grepl(grp_regexp, cross_diff$group) & cross_diff$coefficient == "p",
"value"] <- SEs[SEs[, "parameter"] == param, "p"]
} else { ### if reference level is a subgroup
# create variable to select right rows in jk2 output
#param_finder <- cross_diff$group
#param_finder <- sapply(strsplit(param_finder, ".vs."), function(x) x[1])
## complicated to find param match, because the factor string of another level can contain the string of the current level!
param_selector <- paste0("^", esc_param, "\\.|", "^", esc_param, "_|", "_", esc_param, "\\.|", "_", esc_param, "_")
col_names <- SE_correction[[i]]$refGrp[, "groupName"]
col_levels <- SE_correction[[i]]$refGrp[, "groupValue"]
#if(identical(col_levels, c("female", "TRUE"))) browser()
# filter relevant rows for flexibel number of filter variables
filt_var <- recursive_filter(df = cross_diff, vars = col_names, var_levels = col_levels)
compare_point_estimates(old_est = cross_diff[which(filt_var & grepl(param_selector, cross_diff$group) & cross_diff$coefficient == "est"), "value"],
new_est = SEs[SEs[, "parameter"] == param, "est"],
param = param)
cross_diff[which(filt_var & grepl(param_selector, cross_diff$group) & cross_diff$coefficient == "se"),
"value"] <- SEs[SEs[, "parameter"] == param, "se"]
cross_diff[which(filt_var & grepl(param_selector, cross_diff$group) & cross_diff$coefficient == "p"),
"value"] <- SEs[SEs[, "parameter"] == param, "p"]
}
}
}
# change modus in reporting output for WEC
cross_diff$modus <- unique(SE_correction[[1]][[1]][[1]]$modus)
## mehr Fragen
# _ in Faktorlevels erlaubt?
# next: trend, groupdiffs checken, ob hier SEs korrekt beruecksichtigt werden
# was passiert eig. wenn faktorlevel verschiedener Variablen gleich heissen?
rbind(sep_jk2[["no_cross_diff"]], cross_diff)
}
## Prepare jk2 object for se_correction
separate_jk2 <- function(jk2) {
## Trend or no trend
year <- as.character(unique(jk2[["year"]]))
if(length(year) == 0) year <- 1
## dissect results so far
no_cross_diff <- jk2[is.na(jk2$comparison) | jk2$comparison != "crossDiff" | jk2$parameter != "mean", ]
cross_diff <- jk2[which(jk2$comparison == "crossDiff" & jk2$parameter == "mean"), ]
stopifnot(identical(nrow(no_cross_diff) + nrow(cross_diff), nrow(jk2)))
stopifnot(nrow(cross_diff) > 0)
list(cross_diff = cross_diff, no_cross_diff = no_cross_diff, year = year)
}
recursive_filter <- function(df, vars, var_levels) {
stopifnot(length(vars) == length(var_levels))
filt_var <- rep(TRUE, nrow(df))
for(i in seq_along(vars)) {
filt_var <- filt_var & df[[vars[i]]] == var_levels[[i]]
}
# output: vector with T/F
filt_var
}
seCorrect.rep_se_correction <- function( SE_correction, jk2, grpv ) {
sep_jk2 <- separate_jk2(jk2 = jk2)
cross_diff <- sep_jk2[["cross_diff"]]
no_cross_diff <- sep_jk2[["no_cross_diff"]]
for(i in seq_along(SE_correction)) {
output <- SE_correction[[i]][["resT"]][[sep_jk2[["year"]]]]
rows <- length(SE_correction[[i]][["vgl"]][["groups.divided.by"]])
single_grpv <- SE_correction[[i]][["focGrp"]]
for(param in output[["group"]]) {
compare_point_estimates(old_est = cross_diff[cross_diff$group == param & cross_diff$coefficient == "est", "value"],
new_est = output[output$group == param & output$coefficient == "est", "value"],
param = param)
cross_diff[cross_diff$group == param & cross_diff$coefficient == "se", "value"] <- output[output$group == param & output$coefficient == "se", "value"]
cross_diff[cross_diff$group == param & cross_diff$coefficient == "p", "value"] <- output[output$group == param & output$coefficient == "p", "value"]
}
}
rbind(no_cross_diff, cross_diff)
}
compare_point_estimates <- function(old_est, new_est, param) {
if(abs(abs(old_est) - abs(new_est)) >= 1e-10) {
#browser()
warning("Difference in point estimate of cross level difference for comparison ", param, ": ", round(abs(old_est) - abs(new_est), digits = 3))
}
return()
}
### Hier kann man nicht subset() nehmen, weil es sonst Warnungen beim Paketebauen gibt,
### da subset() die Spaltennamen unquoted haben will
reduceDoubleN <- function(jk2){
cases <- unique(jk2[which(jk2[,"parameter"] == "Ncases"),]) ### cases <- unique(subset(jk2, parameter=="Ncases"))
jk2 <- rbind(jk2[which(jk2[,"parameter"] != "Ncases"),], cases) ### jk2 <- rbind(subset(jk2,parameter != "Ncases"), cases)
return(jk2)}
#reshapeReport <- function(inp, tv, fun, repFunOut, target, wholeGroupName) {
### testweise initialisieren
# inpR <- inp ### untere zeile: levels der Trendvariablen
# lev <- paste(as.character(unique(do.call("rbind", repFunOut[["resT"]])[,repFunOut[["allNam"]][["trend"]]])), collapse="|")
### reshapen nur machen, wenn Parameterspalte mehr als einen Eintrag hat ... das ist bspw. nicht der Fall, wenn man repTable() mit dichotomer AV (Regelstandard erreicht) macht
# if (length(unique(inp[,"parameter"])) > 1) {
### fuer trend
# if (!is.null(tv) ) {
### fuer mittelwerte mit trend
# if (fun == "mean") {
# inpR <- data.frame(tidyr::pivot_wider(inp, names_from = "parameter", values_from = grep(lev, colnames(inp), value=TRUE)), stringsAsFactors=FALSE)
# }
### fuer haeufigkeitsverteilungen mit trend
# if (fun == "table") {
# inpR <- data.frame(tidyr::pivot_wider(inp, names_from = "parameter", values_from = grep(lev, colnames(inp), value=TRUE)), stringsAsFactors=FALSE)
# }
# } else {
### fuer mittelwerte ohne trend
# if (fun == "mean") {
# inpR <- data.frame(tidyr::pivot_wider(inp, names_from = "parameter", values_from = intersect(c("es", "est", "p", "se"), colnames(inp))), stringsAsFactors=FALSE)
# }
### fuer haeufigkeitsverteilungen ohne trend
# if (fun == "table") {
# inpR <- data.frame(tidyr::pivot_wider(inp, names_from = "parameter", values_from = intersect(c("est", "p", "se"), colnames(inp))), stringsAsFactors=FALSE)
# }
# }
# }
### fuer allen Output: spalten, die nur NA haben, loeschen
# inpR <- inpR[, colSums(is.na(inpR)) != nrow(inpR)]
### crossDiff-Benennung anpassen (Mail Juan Jose, 27.04.2022, 9.37 Uhr). Das passiert auch dann, wenn nicht gereshapet wird
# if ("comparison" %in% colnames(inpR)) {
# if( any(c("crossDiff", "trendDiff_cross") %in% inpR[,"comparison"])) {
# rows <- intersect(intersect(grep("wholeGroup$", inpR[,"group"]), setdiff(1:nrow(inpR), grep("_",inpR[,"group"]))), which(inpR[,"comparison"] %in% c("crossDiff", "trendDiff_cross")))
# if (length(rows)>0) {
# inpR[rows,"group"] <- eatTools::removePattern(inpR[rows,"group"], "\\.vs\\.wholeGroup$")
# }
# }
# }
### Alternativ: Mail Kristoph, 12.05.2022, 18.40 Uhr (geht wahrscheinlich nur fuer eine Gruppierungsvariable)
# if (fun == "table" && length(unique(inpR[,"group"])) < nrow(inpR) && target == "BT2021.table" && length(repFunOut[["allNam"]][["group"]])>0 ) {
# prz <- grep("^est|^se", colnames(inpR), value=TRUE)
# if(length(prz)<1) {warning("Cannot found colums with percentage values. Skip value transformation.")} else {
# for ( i in prz) { inpR[,i] <- 100*inpR[,i]}
# }
### Spalten identifizieren, die fuers reshapen entfernt werden muessen, und zwar sowohl fuer mit als auch fuer ohne trend
# weg <- grep(repFunOut[["allNam"]][["group"]], colnames(inpR))
# stopifnot(length(weg)>0)
### erstmal ohne Trend
# if (is.null(tv) ) {
# inpR <- data.frame(tidyr::pivot_wider(inpR[,-weg], names_from = "comparison", values_from = intersect(c("es", "est", "p", "se"), colnames(inpR))), stringsAsFactors=FALSE)
# } else { ### mit Trend
# inpR <- data.frame(tidyr::pivot_wider(inpR[,-weg], names_from = "comparison", values_from = grep(lev, colnames(inp), value=TRUE)), stringsAsFactors=FALSE)
# }
# inpR[,"group"] <- car::recode(inpR[,"group"], paste0("'wholeGroup'='",wholeGroupName,"'"))
# colnames(inpR) <- eatTools::removePattern(colnames(inpR), pattern = "_NA$")
# }
# inpR <- inpR[, colSums(is.na(inpR)) != nrow(inpR)] ### nochmal spalten loeschen mit nur NA
# return(inpR)}
Try the eatRep package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.