R/bmdFunctions_visualization.R
In alfcrisci/bmdModeling: Benchmark dose modeling
#' Wrapper function for plotting the results of proast
#' @param ans.all list, with all results that were obtained during the analysis
#' @param track logical, if TRUE (FALSE by default) print the name of the function which is currently being run
#' @param plot.type, integer, defines the type of plot; for cont choose 1 (y vs x),
#' for dtype 2 or 3 choose one of 3 (y.up vs. x), 5 (cumulative y vs. x);
#' default is 1 for cont and 5 else
#' @return list, updated version of ans.all
#' @export
f.plot.all <- function (ans.all, track = FALSE, plot.type = NA) {
if (track) print("f.plot.all")
# MV added default plot type choice
if(is.na(plot.type)){
if(ans.all$cont){
ans.all$plot.type <- 1
} else {
ans.all$plot.type <- 5
}
} else {
ans.all$plot.type <- plot.type
}
if (is.na(ans.all$xy.lim[1]) || ans.all$xy.lim[1] == 0) {
x <- ans.all$x
dum.contr <- min(x[x > 0])/5
ans.all$xy.lim <- c(dum.contr, min(x), max(x), min(ans.all$y), max(ans.all$y))
}
with(ans.all, {
if (cont) {
ans.all$y.lim <- f.plot.con(ans.all, track = track)
if (model.ans != 11){
f.lines.con(ans.all, track = track)
}
if (!is.na(CED[1]) & model.ans != 16) {
f.cedlines.con(ans.all, track = track)
}
} else {
if (dtype %in% c(4, 6, 84) | (dtype == 2 & plot.type < 3)) {
ans.all <- f.plot.frq(ans.all, track = track)
if (model.ans != 14)
ans.all$gr.txt <- f.lines.frq(ans.all, track = track)
f.cedlines.bin(ans.all, track = track)
} else if(dtype %in% 2:3) {
## in f.mm4.cat, code without combi:
ans.all$combi <- F
ans.all$categ.ans <- 0
# f.plot.all(ans.all)
# ans.all.plt <- ans.all #f.model.bb(ans.all = ans.all)
# Error when plotting lines. In Proast61.3, NOTE: Full model may not be applicable for dtype = 2
if (!(model.ans == 14 & dtype == 2))
ans.all$y.lim <- f.lines.cat(ans.all, track = track)
title(ans.all$model.name)#f.mtext(ans.all)
# in f.cat: code with combi.ans:
## ans.all$combi.ans <- 1
## ans.all$categ.ans <- 0
## ans.all$plot.type <- 5
## f.plot.all(ans.all, sep = F)
## f.graph.window(4) equivalent to:
# par(mfcol = c(2, 2), mar = c(3.2, 3.7, 2, 1.5), cex.main = 1.2,
# cex.sub = 0.9, cex.lab = 0.9, cex = 0.9)
# ans.all.plt <- ans.all
# ans.all.plt$modelname <- ""
# ans.all.plt$plot.type <- 3
# ans.all.plt$ced.lines <- TRUE
# ans.all.plt$main <- ""
# ans.all.plt$y.lim <- f.lines.cat(ans.all.plt, track = track)
# title(model.name, col = color[1], cex.main = 0.9, font.main = 1)
## f.mtext(ans.all.plt, track = track)
# ans.all.plt$plot.type <- 9
# ans.all.plt$modelname <- ""
# ans.all.plt$ced.lines <- F
# ans.all.plt$y.lim <- f.lines.cat(ans.all.plt, track = track)
#
# if (!is.na(CED.matr[1])) {
# th.par.cum <- cumsum(th.par)
# CED.vect <- CED.matr[1, ]
# if (ans.scale == 1) {
# CED.vect <- log10(CED.vect)
# x <- log10(x[x > 0]/sf.x)
# }
# for (ii in 1:length(CED.vect))
# lines(c(min(x), CED.vect[ii]), rep(th.par.cum[ii], 2), lty = 2)
# y.leg.2 <- paste(y.leg, " (latent var.uu)")
# if (dtype == 3 & nr.gr == 1) {
# CED.tmp <- round(CED.vect, rep(4, length(CED.vect)))
# CED.tmp <- CED.tmp[CED.tmp < max(x)]
# CED.left <- c(CED.tmp[1]/100, CED.tmp)
# CED.right <- c(CED.tmp, max(x))
# CED.mid <- (CED.left + CED.right)/2
# if (ans.scale == 1)
# CED.mid <- 10^(CED.mid)
# uu <- f.expect.con(model.ans, CED.mid, regr.par,
# fct1 = 1, fct2 = 1, CES = CES, track = track)
# zz <- logb(uu)
# if (ans.scale == 1)
# CED.mid <- log10(CED.mid)
# for (ii in 1:length(scores.orig))
# text(CED.mid[ii], zz[ii], scores.orig[ii])
# }
# }
# th.par.cum <- cumsum(th.par)
# if (mode(CED.matr) == "NULL")
# CED.vect <- 0
# else if (is.na(CED.matr[1]))
# CED.vect <- 0
# else CED.vect <- CED.matr[1, ]
# ans.all.plt$plot.type <- 5
}
}
if (track)
print("f.plot.all : END")
return(ans.all)
})
}
#' Plot the predicted response values (points) and estimated CIs
#' @param ans.all list, with all results that were obtained during the analysis
#' @param track logical, if TRUE (FALSE by default) print the name of the function which is currently being run
#' @return vector, range of the CI bounds on the predicted response values
#' @export
f.plot.con <- function (ans.all, track = FALSE) {
if (track)
print("f.plot.con")
with(ans.all, {
if (is.null(x.leg))
x.leg <- ""
if (is.null(y.leg))
y.leg <- ""
nr.aa <- max(fct1)
nr.bb <- max(fct2)
nr.var <- max(fct3)
twice <- FALSE
if (sum(fct1 != fct2) == 0)
twice <- TRUE
f.means <- function(xx, yy, dtype, sd2.log, plot.type, nn) {
if (track)
print("f.means within f.plot.con")
mean.y <- yy
if (dtype %in% c(1, 5, 15, 101)){
mean.y <- exp(tapply(log(yy), xx, mean))
} else if (dtype == 26){
mean.y <- tapply(sqrt(yy), xx, mean)^2
} else if (dtype == 25){
mean.y <- tapply(yy, xx, mean)
}
if (plot.type %in% c(1, 2, 5, 10, 11))
mean.y <- mean.y
if (plot.type %in% c(3, 4, 6))
mean.y <- log10(mean.y)
if (plot.type %in% c(7, 8))
mean.y <- sqrt(mean.y)
f.var <- function(dtype, xx, yy, sd2.log, nn) {
if (track)
print("f.var within f.plot.con")
if (dtype %in% c(10, 15, 110, 250, 260)) {
SS <- sum(sd2.log * (nn - 1))
var.within <- SS/sum(nn - 1)
df <- sum(nn - 1)
} else {
if (dtype %in% c(1, 5)) {
yy <- log(yy)
} else if (dtype == 26){
yy <- sqrt(yy)
}
mn <- tapply(yy, xx, mean)
nn <- tapply(yy, xx, length)
Vmn <- rep(mn[1], nn[1])
if (dtype %in% c(1, 5)){
for (ii in 2:length(mn)){
Vmn <- c(Vmn, rep(mn[ii], nn[ii]))
}
}
resid <- yy - Vmn
df <- sum(nn - 1)
var.within <- (sum(nn) - 1) * var(resid)/df
}
out.lst <- list(Vsem = sqrt(var.within/nn), df = df)
return(out.lst)
}
var.out <- f.var(dtype, xx, yy, sd2.log, nn)
Vsem <- var.out$Vsem
df <- var.out$df
Vconf <- qt(0.975, df) * Vsem
if (dtype %in% c(1, 5, 10, 15)) {
conf.L <- mean.y/exp(Vconf)
conf.U <- mean.y * exp(Vconf)
} else if (dtype %in% c(25, 250)) {
conf.L <- mean.y - Vconf
conf.U <- mean.y + Vconf
} else {
conf.L <- (sqrt(mean.y) - Vconf)^2
conf.U <- (sqrt(mean.y) + Vconf)^2
}
y.lim.tmp <- range(c(conf.L, conf.U))
if (plot.type %in% c(3, 4, 6)) {
conf.L <- log10(conf.L)
conf.U <- log10(conf.U)
}
if (plot.type %in% c(7, 8)) {
conf.L <- sqrt(conf.L)
conf.U <- sqrt(conf.U)
}
y.lim.tmp <- c(min(conf.L), max(conf.U))
out.lst <- list(conf.L = conf.L, conf.U = conf.U,
y.lim.tmp = y.lim.tmp, mean.y = mean.y)
return(out.lst)
}
x.plt <- x
y.plt <- yy
switch(as.character(plot.type),
'1' = { #1
x.plt <- x
y.plt <- yy
},
'3' = { #3
x.plt <- x
y.plt <- log10(yy)
},
'5' = { #5
y.plt <- yy
x.plt <- sqrt(x)
if (x.leg != "") xleg <- paste("sqrt-", x.leg, sep = "") else xleg <- ""
})
if (dtype %in% c(10, 110, 250, 260)) {
plt.mns <- 3
} else {
plt.mns <- 1
}
cex.1 <- 0.6
cex.2 <- 1.5
mark <- c(1:2, 4:25, 33:500)
# if (nr.aa == 1 && nr.bb == 1 && nr.var == 1) {
#
# if (max(sp.fact) > 1) {
#
# fct2 <- sp.fact
# nr.bb <- max(fct2)
# cat("\n\nno covariates, but subgroups are plotted distinctly\n\n")
#
# }
#
# }
calculationsPlot <- f.means(xx = x.plt, yy = y.plt, dtype = dtype,
sd2.log = sd2.log, plot.type = plot.type, nn = nn)
yRange <- range(c(y.plt, calculationsPlot$y.lim.tmp))
plot(x.plt, y.plt, main = heading, ylim = yRange,
xlab = x.leg, ylab = y.leg, type = "n", col = color[1])
if (nr.aa == 1 & nr.bb == 1 & nr.var > 1) {
nr.bb <- nr.var
fct2 <- fct3
}
if (nr.aa > 1 & nr.bb > 1 && twice) {
nr.aa <- 1
fct1 <- rep(1, length(x))
}
shift.tmp <- 0
shift.abs <- 0
bbb <- (max(x.plt, na.rm = TRUE) - min(x.plt, na.rm = TRUE))/70
zz <- 0
for (jj in 1:nr.bb) {
for (ii in 1:nr.aa) {
x.part <- x.plt[fct1 == ii & fct2 == jj]
y.part <- y.plt[fct1 == ii & fct2 == jj]
y.tmp <- yy[fct1 == ii & fct2 == jj]
if (length(y.part) > 0) {
zz <- zz + 1
if (zz > 500)
cat("\nATTENTION: number of subgroups too large for plotting\n\n")
type.dum <- "p"
lt.y <- 1
if (model.ans == 11) {
lt.y <- 2
type.dum <- "b"
}
x.part <- x.part + shift.tmp
if (plt.mns < 3) {
points(x.part, y.part, col = color[zz], pch = mark[zz],
cex = cex.1, type = "p", lty = lt.y)
}
mean.x <- NA
# Plot mean
# Note: changed compared to original code, because
# Error in xy.coords(x, y) : 'x' and 'y' lengths differ
if (dtype %in% c(1, 5, 15, 25, 26)) {
mean.x <- as.numeric(tapply(x.part, x.part, mean))
mean.y <- f.means(x.part, y.tmp, dtype, plot.type = plot.type)$mean.y
points(mean.x, mean.y, col = color[zz], pch = mark[zz],
cex = 1, type = type.dum, lty = lt.y)
} else {
mean.x <- x.part
mean.y <- y.part
points(mean.x, mean.y, col = color[zz], pch = mark[zz],
cex = 1, type = type.dum, lty = lt.y)
}
if (dtype %in% c(10, 250, 260)) {
sd2.log.part <- sd2.log[fct1 == ii & fct2 == jj]
nn.part <- nn[fct1 == ii & fct2 == jj]
}
out.lst <- f.means(x.part, y.tmp, dtype, sd2.log = sd2.log.part,
plot.type, nn = nn.part)
conf.L.part <- out.lst$conf.L
conf.U.part <- out.lst$conf.U
# Plot confidence limits
for (qq in 1:length(conf.L.part)) {
mean.x[qq] <- mean.x[qq]
lines(rep(mean.x[qq], 2), c(conf.L.part[qq],
conf.U.part[qq]), col = color[zz])
lines(c(mean.x[qq] - bbb, mean.x[qq] +
bbb), rep(conf.L.part[qq], 2), col = color[zz])
lines(c(mean.x[qq] - bbb, mean.x[qq] +
bbb), rep(conf.U.part[qq], 2), col = color[zz])
}
}
shift.tmp <- shift.tmp + shift.abs
}
}
ans.all$y.lim <- range(y.plt)
ans.all$y.lim.CI <- range(y.plt)
if (track)
print("f.plot.con: END")
return(calculationsPlot$y.lim.tmp)
})
}
#' Plot dotted lines for the determined CED value for continous model
#' @param ans.all list, with all results that were obtained during the analysis
#' @param track logical, if TRUE (FALSE by default) print the name of the function which is currently being run
#' @return NULL
#' @export
f.cedlines.con <- function (ans.all, track = FALSE) {
if (track) print("f.cedlines.con")
with(ans.all, {
nr.subgr <- nrow(regr.par.matr)
dum.contr <- xy.lim[1]
# MV changed this
# low.y <- min(y.lim)
low.y <- xy.lim[4]
for (jj in 1:nr.subgr) {
CED.0 <- unlist(CED[jj])
CES.0 <- CES
ES.0 <- f.expect.con(model.ans, CED.0, regr.par.matr[jj,],
fct1 = 1, fct2 = 1, fct5 = 1, CES = CES.0,
increase = increase, track = track)
ES.y <- rep(ES.0, 2)
ES.x <- c(min(x), CED.0)
CED.y <- c(low.y, ES.0)
CED.x <- rep(CED.0, 2)
lines(ES.x, ES.y, lty = 2, lwd = 1)
lines(CED.x, CED.y, lty = 2, lwd = 1)
}
if (track) print("f.cedlines.con: END")
return(NULL)
})
}
#' Plot the curve of the estimated model
#' @param ans.all list, with all results that were obtained during the analysis
#' @param track logical, if TRUE (FALSE by default) print the name of the function which is currently being run
#' @return NULL
#' @export
f.lines.con <- function (ans.all, track = FALSE) {
if (track)
print("f.lines.con")
with(ans.all, {
f.lines.tmp <- function(ans.all.tmp) {
with(ans.all.tmp, {
nbins <- 1000
nbins.0 <- nbins/10
xline <- seq(from = min(x, na.rm = TRUE),
to = max(x, na.rm = TRUE), length = nbins)
twice <- FALSE
expect <- f.expect.con(model.ans, xline, regr.par,
fct1 = rep(1, length(xline)), fct2 = rep(1,
length(xline)), fct3 = 1, CES = CES, twice = twice,
increase = increase, x.mn = x.mn, track = track)
lines(xline, expect, lty = 1, col = col.tmp)
return(NULL)
})
}
ans.all.tmp <- ans.all
ans.all.tmp$dum.contr <- xy.lim[1]
ans.all.tmp$min.x <- xy.lim[2]
ans.all.tmp$max.x <- xy.lim[3]
if ((nr.aa > 1 || nr.bb > 1) && (nr.cc > 1 || nr.dd > 1) && !twice) {
warning("Combination of covariate on c or d and yet another covariate
\ndoes not allow plotting of all individual curves")
}
if (ans.plt == 1)
regr.par.matr <- matrix(regr.par, nrow = 1)
else regr.par.matr <- f.pars(ans.all.tmp)$regr.par.matr
count <- 0
if (nr.aa == 1 && nr.bb == 1 && nr.var > 1){
for (ii in 1:nr.var) {
count <- count + 1
ans.all.tmp$regr.par <- regr.par.matr[1, ]
ans.all.tmp$col.tmp <- color[count]
f.lines.tmp(ans.all.tmp)
}
} else {
kk.max <- nrow(regr.par.matr)
par.bb.old <- NA
par.bb.new <- 0
for (kk in 1:kk.max) {
ans.all.tmp$regr.par <- regr.par.matr[kk, ]
if (model.ans > 1)
par.bb.new <- ans.all.tmp$regr.par[2]
if (ans.plt == 1) {
if (!identical(par.bb.new, par.bb.old) || (model.ans ==
1)) {
count <- count + 1
ans.all.tmp$col.tmp <- color[count]
}
}
else {
count <- count + 1
ans.all.tmp$col.tmp <- color[count]
}
f.lines.tmp(ans.all.tmp)
par.bb.old <- par.bb.new
}
}
if (track)
print("f.lines.con : END")
return(NULL)
})
}
#' Plot the predicted response values (points) and estimated CIs
#' @param ans.all list, with all results that were obtained during the analysis
#' @param track logical, if TRUE (FALSE by default) print the name of the function which is currently being run
#' @return updated ans.all
#' @export
f.plot.frq <- function (ans.all, track = FALSE){
if(track) print("f.plot.frq")
with(ans.all, {
if (length(fct1) == 1)
fct1 <- rep(1, length(x))
if (length(fct2) == 1)
fct2 <- rep(1, length(x))
dum.contr <- xy.lim[1]
y.lim <- c(xy.lim[4], xy.lim[5])
x.lim <- c(xy.lim[2], xy.lim[3])
if (down) y <- 1 - y
if (max(sp.fact) > 1) fct2 <- sp.fact
if (max(fct1) > 1 & max(fct2) > 1 & (sum(fct1 == fct2) == length(fct1)))
twice <- T
shift.abs <- 0
shift.tmp <- 0
mark <- c(1:25, 32:255)
x.lim.plt <- x.lim
y.lim.plt <- y.lim
if (shift > 0) shift.abs <- (max(x) - min(x))/shift
fct1 <- as.numeric(fct1)
fct2 <- as.numeric(fct2)
dum <- max(fct1, fct2)
plot(x, y, #ylim = y.lim, xlim = x.lim, #removed because error if NA
main = heading, xlab = x.leg, ylab = y.leg, type = "n",
col = color[1])
if (shift == 0)
bbb <- (max(x) - min(x))/70 else
bbb <- (max(x) - min(x))/shift
if (max(as.numeric(fct3)) > 1) fct1 <- fct3
zz <- 0
for (jj in 1:max(fct2)) {
for (ii in 1:max(fct1)) {
x.tmp <- x[fct1 == ii & fct2 == jj]
if (!is.na(x.tmp[1])) {
zz <- zz + 1
x.tmp <- x.tmp + shift.tmp
y.tmp <- y[fct1 == ii & fct2 == jj]
n.tmp <- nn[fct1 == ii & fct2 == jj]
if (is.na(cex.1)) {
if (dtype == 6)
cex.1 <- 0.6
else cex.1 <- 1.2
}
points(x.tmp, y.tmp, col = color[zz], pch = mark[zz], cex = cex.1)
if (dtype == 4 & model.ans == 14)
lines(x.tmp, y.tmp, col = color[zz], lty = 2)
if (CI.plt) {
k.tmp <- kk[fct1 == ii & fct2 == jj]
k.tmp <- k.tmp[!is.na(k.tmp)]
n.tmp <- n.tmp[!is.na(n.tmp)]
L025 <- numeric()
L975 <- numeric()
for (qq in 1:length(k.tmp)) {
LL.tmp <- f.LL.bin(k.tmp[qq], n.tmp[qq])
L025[qq] <- LL.tmp[1]
L975[qq] <- LL.tmp[2]
}
if (dtype %in% c(4, 6, 84)) {
points(x.tmp, L025, pch = "-", cex = 1.2, col = color[zz])
points(x.tmp, L975, pch = "-", cex = 1.2, col = color[zz])
for (qq in 1:length(L025)) {
lines(rep(x.tmp[qq], 2), c(L025[qq], L975[qq]), col = color[zz])
lines(c(x.tmp[qq] - bbb, x.tmp[qq] + bbb), rep(L025[qq], 2), col = color[qq])
lines(c(x.tmp[qq] - bbb, x.tmp[qq] + bbb), rep(L975[qq], 2), col = color[zz])
}
}
shift.tmp <- shift.tmp + shift.abs
}
}
}
}
if (!is.na(fct1.full[1]))
if (dtype == 6) {
if (is.null(ans.all$cex.2))
cex.2 <- 1.5
zz <- 0
for (jj in 1:max(fct2.full)) {
for (ii in 1:max(fct1.full)) {
x.tmp <- x.full[fct1.full == ii & fct2.full ==
jj]
if (!is.na(x.tmp[1])) {
zz <- zz + 1
x.tmp <- x.tmp + shift.tmp
y.tmp <- pi.full[fct1.full == ii & fct2.full ==
jj]
if (sum(is.na(pi.full)) == 0) {
if (model.ans == 14)
type.dum <- "b"
else type.dum <- "p"
points(x.tmp, y.tmp, cex = cex.2, pch = mark[zz],
col = color[zz], type = type.dum, lty = 2)
}
}
}
}
}
if (track) print("plot.frq: END")
ans.all$y.lim <- y.lim
ans.all$y.lim.plt <- y.lim.plt
ans.all$shift.abs <- shift.abs
return(ans.all)
})
}
#' Plot the lines
#' @param ans.all list, with all results that were obtained during the analysis
#' @param track logical, if TRUE (FALSE by default) print the name of the function which is currently being run
#' @return ans.all$gr.txt
#' @export
f.lines.frq <- function (ans.all, track = FALSE){
if (track) print("f.lines.frq")
with(ans.all, {
shift.abs <- 0
if (model.ans == 14) return(invisible())
x.high <- max(x)
f.lines.tmp <- function(model.ans, x, regr.par,
th.par, sig.par, plot.type, dtype, l.ty = l.ty,
color, x.high, shift, CES = CES,
twice, pi.full, down) {
if (track) print("f.lines.tmp in f.lines.frq")
sat.mod <- FALSE
if (model.ans == 14) sat.mod <- TRUE
if (!is.na(pi.full[1]) & sat.mod) {
if (dtype %in% c(4, 84))
xline <- x
f(dtype == 6)
xline <- tapply(x, x, mean)
expect <- pi.full
expect <- expect[order(xline)]
xline <- sort(xline)
} else {
xline <- seq(from = min(x), to = x.high, length = 1000)
expect <- f.expect.bin(model.ans, x = xline, regr.par,
CES = CES, ttt = ttt, twice = twice, ces.ans = ces.ans,
CES1 = CES1, CES2 = CES2)
expectdum <- f.expect.bin(model.ans, x = 0, regr.par,
CES = CES, ttt = ttt, twice = twice, ces.ans = ces.ans,
CES1 = CES1, CES2 = CES2)
if (down) {
expect <- 1 - expect
expectdum <- 1 - expectdum
}
pipi <- f.expect.bin(model.ans, x = x, regr.par,
CES = CES, ttt = ttt, twice = twice, ces.ans = ces.ans,
CES1 = CES1, CES2 = CES2)
if (sum(is.na(expect)) > 0) return()
xline <- xline + shift
}
lines(xline, expect, lty = l.ty, col = color)
}
if (max(as.numeric(fct3)) > 1) {
ans.all <- f.pars.frq(ans.all, track = track)
regr.par.matr <- ans.all$regr.par.matr
th.par.vec <- ans.all$th.par.vec
kk.max <- length(th.par.vec)
# MV added this to prevent errors for model.ans = 1 with covariates
if(kk.max != nrow(regr.par.matr))
regr.par.matr <- t(regr.par.matr)
} else {
ans.all <- f.pars(ans.all, track = track)
regr.par.matr <- ans.all$regr.par.matr
# MV changed this to prevent showing too many lines
# kk.max <- length(regr.par.matr[, 1])
kk.max <- nrow(unique(regr.par.matr))
th.par.vec <- rep(th.par, kk.max)
}
shift.tmp <- 0
for (kk in 1:kk.max) {
f.lines.tmp(model.ans, x, regr.par.matr[kk, ],
th.par = th.par.vec[kk], sig.par, plot.type,
dtype = dtype, l.ty = l.ty,
col = color[kk], x.high = x.high,
shift = shift.tmp, CES = CES,
twice = twice, pi.full = pi.full,
down = down)
shift.tmp <- shift.tmp + shift.abs
}
if (track) print("f.lines.frq : END")
return(ans.all$gr.txt)
})
}
#' ?
#' @param ans.all list, with all results that were obtained during the analysis
#' @param track logical, if TRUE (FALSE by default) print the name of the function which is currently being run
#' @return updated ans.all
#' @export
f.pars.frq <- function (ans.all, track = FALSE)
{
if (track) print("f.pars.frq")
with(ans.all, {
if (max(fct1) == 1 & max(fct2) == 1 & max(as.numeric(fct3)) == 1)
regr.par.matr <- matrix(regr.par, nrow = 1)
else {
if (max(as.numeric(fct3)) == 1)
pars <- c(regr.par, th.par)
else {
th.lvm <- f.th.lvm(th.par, fct3, fct3.ref, track = track)
pars <- c(regr.par, th.lvm)
fct3 <- as.numeric(factor(fct3))
}
nrp <- length(regr.par)
regr.par.matr <- numeric()
nr.aa <- max(fct1)
nr.bb <- max(fct2)
nr.th <- max(as.numeric(fct3))
if (!covar.dd) {
if ((model.ans == 11) & (nr.th > 1))
fct1 <- fct3
par.tmp <- rep(NA, length(pars))
if (nrp > (nr.aa + nr.bb))
par.rest <- pars[(nr.aa + nr.bb + 1):nrp]
else par.rest <- numeric()
zz <- 1
for (jj in 1:nr.bb) {
f1 <- fct1[fct2 == jj]
f1.lev <- levels(factor(f1))
for (ii.index in f1.lev) {
ii <- as.numeric(ii.index)
if (max(fct3) > 1) {
f3 <- fct3[fct1 == ii & fct2 == jj]
f3.lev <- levels(factor(f3))
for (kk.index in f3.lev) {
kk <- as.numeric(kk.index)
if (!twice)
gr.txt[zz] <- paste(fct2.txt[jj], fct3.txt[kk], sep = "-")
if (twice)
gr.txt[zz] <- fct2.txt[jj]
par.tmp <- c(pars[ii], pars[nr.aa + jj],
par.rest, pars[nr.aa + nr.bb + length(par.rest) + kk])
regr.par.matr <- rbind(regr.par.matr, par.tmp)
zz <- zz + 1
}
}
if (max(fct3) == 1) {
par.tmp <- c(pars[ii], pars[nr.aa + jj], par.rest, th.par)
regr.par.matr <- rbind(regr.par.matr, par.tmp)
if (!twice)
gr.txt[kk] <- paste(fct1.txt[ii], fct2.txt[jj],
sep = "-")
if (twice)
gr.txt[kk] <- fct1.txt[ii]
}
}
}
n.col <- nrp - nr.aa - nr.bb + 2 + 1
regr.par.matr <- matrix(regr.par.matr[, 1:n.col],
ncol = n.col)
}
if (covar.dd) {
cc <- pars[nr.aa + nr.bb + 1]
for (jj in 1:nr.bb) {
par.tmp <- c(pars[jj], pars[nr.aa + jj], cc,
pars[nr.aa + nr.bb + 1 + jj])
regr.par.matr <- rbind(regr.par.matr, par.tmp)
}
}
}
npar <- length(regr.par.matr[1, ])
th.par.vec <- regr.par.matr[, npar]
regr.par.matr <- regr.par.matr[, 1:(npar - 1)]
# MV added this to prevent errors
if(!is.matrix(regr.par.matr))
regr.par.matr <- matrix(regr.par.matr, nrow = 1)
ans.all$regr.par.matr <- regr.par.matr
ans.all$gr.txt <- gr.txt
ans.all$th.par.vec <- th.par.vec
if (track) print("f.pars.frq: END")
return(ans.all)
})
}
#' compute th.lvm?
#' @param th.par parameter theta
#' @param fct3 initial value for theta
#' @param fct3.ref theta reference?
#' @param track logical, if TRUE (FALSE by default) print the name of the function which is currently being run
#' @return th.lvm
#' @export
f.th.lvm <- function (th.par, fct3, fct3.ref, track = FALSE) {
if (track) print("f.th.lvm")
nr.th <- length(levels(as.factor(fct3)))
fct3.lev <- as.numeric(levels(as.factor(fct3)))
th.lvm <- rep(th.par, nr.th)
th.lvm <- th.lvm - (log(fct3.ref) - log(fct3.lev))
if (track) print("f.th.lvm: END")
return(th.lvm)
}
#' Plot dotted lines for the determined CED value for categorical model
#' @param ans.all list, with all results that were obtained during the analysis
#' @param inputResponse numeric, value for the response at estimated CED;
#' default value is NA
#' @param inputCED numeric, value for the CED; default value is NA
#' @param track logical, if TRUE (FALSE by default) print the name of the function which is currently being run
#' @return NULL
#' @export
f.cedlines.bin <- function (ans.all, inputResponse = NA, inputCED = NA,
track = FALSE) {
if (track) print("f.cedlines.bin")
with(ans.all, {
# MV added these conditions to use plot for model averaging
if(is.na(inputResponse))
response <- response.matr else
response <- inputResponse
if(is.na(inputCED))
CED <- CED.matr else
CED <- inputCED
dum.contr <- xy.lim[1]
nr.lev <- max(max(fct1, fct2))
ES.y <- matrix(response[1], ncol = 2)
ES.x <- matrix(c(x[1], CED[1]), ncol = 2)
if (nr.lev > 1)
for (ii in (2:nr.lev)) {
ES.y <- rbind(ES.y, rep(response[ii], 2))
ES.x <- rbind(ES.x, c(x[1], CED[ii]))
}
CED.y <- matrix(c(min(y), response[1]), ncol = 2)
CED.x <- matrix(rep(CED[1], 2), ncol = 2)
if (nr.lev > 1)
for (ii in (2:nr.lev)) {
CED.y <- rbind(CED.y, c(min(y), response[ii]))
CED.x <- rbind(CED.x, rep(CED[ii], 2))
}
for (ii in (1:nr.lev)) {
if (down) {
ES.y[ii, ] <- 1 - ES.y[ii, ]
CED.y[ii, 2] <- 1 - CED.y[ii, 2]
}
lines(ES.x[ii, ], ES.y[ii, ], lty = 2)
lines(CED.x[ii, ], CED.y[ii, ], lty = 2)
}
return(NULL)
})
}
#' Calculate confidence intervals for each of the plotted points for quantal response
#' @param k integer, number of events in given group
#' @param n integer, number of observations in given group
#' @param n.pi ?
#' @param conf ?
#' @param step ?
#' @return vector with lower and upper limit
#' @export
f.LL.bin <- function (k = 4, n = 20, n.pi = 100, conf = 0.9, step = 1.2)
{
compl <- 0
alfa <- (1 - conf)/2
if ((k/n) > 0.5) {
compl <- T
k <- n - k
}
pi.min <- (k/n)
pi.max <- 1
pi.range <- 0.999 * seq(pi.min, pi.max, (pi.max - pi.min)/n.pi)
prob1 <- pbinom(k, n, pi.range)
if (sum(pi.range < 0 | pi.range > 1) > 0) {
}
dum1 <- pi.range[prob1 <= alfa]
while (length(dum1) > n.pi/2) {
pi.min <- pi.min/step
pi.max <- pi.max/step
pi.range <- 0.999 * seq(pi.min, pi.max, (pi.max - pi.min)/n.pi)
prob1 <- pbinom(k, n, pi.range)
if (sum(pi.range < 0 | pi.range > 1) > 0) {
}
dum1 <- pi.range[prob1 <= alfa]
}
dum2 <- pi.range[prob1 >= alfa]
P2 <- (max(prob1[prob1 <= alfa]))
P1 <- (min(prob1[prob1 >= alfa]))
L.up <- ((P1 - alfa) * (min(dum1) - max(dum2)))/(P1 - P2) +
max(dum2)
pi.min <- (k/n)
pi.max <- 1
pi.range <- 0.999 * seq(pi.min, pi.max, (pi.max - pi.min)/n.pi)
if (k == 0)
L.low <- 0
else {
prob2 <- 1 - pbinom(k - 1, n, pi.range)
if (sum(pi.range < 0 | pi.range > 1) > 0) {
}
if (sum(pi.range < 0 | pi.range > 1) > 0) {
}
dum1 <- pi.range[prob2 <= alfa]
while (length(dum1) < n.pi/2) {
pi.max <- pi.max/step
pi.min <- pi.min/step
pi.range <- 0.999 * seq(pi.min, pi.max, (pi.max -
pi.min)/n.pi)
prob2 <- 1 - pbinom((k - 1), n, pi.range)
if (sum(pi.range < 0 | pi.range > 1) > 0) {
}
dum1 <- pi.range[prob2 <= alfa]
}
dum2 <- pi.range[prob2 >= alfa]
P2 <- (max(prob2[prob2 <= alfa]))
P1 <- (min(prob2[prob2 >= alfa]))
L.low <- min(dum2) - ((P1 - alfa) * (min(dum2) - max(dum1)))/(P1 -
P2)
}
if (compl) {
L.low.tmp <- L.low
L.low <- 1 - L.up
L.up <- 1 - L.low.tmp
}
return(c(L.low, L.up))
}
#' plot lines for categorical model
#'
#' Note: used in the f.plot.all function, for dtype == 2, 3,
#' and where plot.type is set to 3, 5, 9?
#' so only this part of the code is retained
#' @param ans.all list, with all results that were obtained during the analysis
#' @param track logical, if TRUE (FALSE by default) print the name of the function which is currently being run
#' @return c(0, y.max)
#' @export
f.lines.cat <- function (ans.all, track = FALSE){
if (track) print("f.lines.cat")
if (length(ans.all$cex.tmp) == 0) cex.tmp <- 1
with(ans.all, {
dum.contr <- xy.lim[1]
x.lim <- c(xy.lim[2], xy.lim[3])
y.max <- xy.lim[5]
f.ced.lines <- function(regr.par, nth, dum.contr, CED,
x, response, plot.type, upper.uu = upper.uu, col.fct) {
if (track) print("f.ced.lines in f.lines.cat")
switch(as.character(plot.type),
'3' = for (k in 1:nth)
lines(rep(CED[k], 2), c(0, response[k]),
lty = 2, col = col.fct),
'9' = for (k in 1:nth)
lines(rep(CED[k], 2), c(min(logb(regr.par[2])), th.par[k]),
lty = 2, col = col.fct)
)
switch(as.character(plot.type),
# 3
'3' = for (k in 1:nth)
lines(c(0, CED[k]), c(response[k], response[k]),
lty = 2, col = col.fct),
# 9
'9' = for (k in 1:nth)
lines(c(0, CED[k]), rep(th.par[k], 2), lty = 2, col = col.fct)
)
# nothing is done for plot.type == 5
}
f.lines.tmp <- function(x, y, regr.par, th.par, sig.par,
nth, dum.contr, plot.type, col.fct, categ.ans = categ.ans,
first.plot, model.ans, model.type, x.leg, y.leg,
x.high, upper.uu = upper.uu, CES.cat, CES, reverse = reverse,
def.exc, cex.tmp, show, dtype, y.max, x.lim, twice) {
cat.txt <- c("category 1", "category 2", "category 3",
"category 4", "category 5", "category 6", "category 7",
"category 8")
if (track) print("f.lines.tmp in f.lines.cat")
xline <- seq(from = dum.contr, to = x.high, length = 100)
if (model.type == 1) {
#MV error when plot for dtype = 2, so added first part using dtype = 2 instead of default dtype = 6
if(dtype == 2){
pipi <- matrix(f.expect.bin(model.ans, x, regr.par,
CES = CES, ces.ans = ces.ans, dtype = dtype, track = track), ncol = 1)
piline <- matrix(f.expect.bin(model.ans, xline,
regr.par, CES = CES, ces.ans = ces.ans, dtype = dtype, track = track), ncol = 1)
} else {
pipi <- matrix(f.expect.bin(model.ans, x, regr.par,
CES = CES, ces.ans = ces.ans, track = track), ncol = 1)
piline <- matrix(f.expect.bin(model.ans, xline,
regr.par, CES = CES, ces.ans = ces.ans, track = track), ncol = 1)
}
} else {
pipi <- f.expect.cat(model.type, model.ans, x,
regr.par, th.par, sig.par, CES = CES, CES.cat = CES.cat,
ces.ans = ces.ans, dtype = dtype, twice = twice, track = track)
piline <- f.expect.cat(model.type, model.ans,
xline, regr.par, th.par, sig.par, CES = CES,
CES.cat = CES.cat, ces.ans = ces.ans, dtype = dtype,
twice = twice, track = track)
}
switch(as.character(plot.type),
# 3
'3' = {
pi.up <- matrix(0, length(xline), (nth + 1))
if (first.plot == T)
plot(xline, pi.up[, 1],
xlim = x.lim, ylim = c(0, 1), xlab = x.leg,
ylab = paste(y.leg, "(fraction of scores)"),
type = "n", col = col.fct)
for (k in 1:nth) {
dum <- nth + 1 - k
pi.up[, dum] <- piline[, dum] + pi.up[, (dum + 1)]
}
for (k in 1:nth)
lines(xline, pi.up[, k], col = col.fct, lty = k)
y.max <- 1
x.max <- x.high
first.plot <- F
},
'5' = {
y.leg <- paste(y.leg, " (cumu)")
# if (!is.na(def.exc[1])) {
# de.lev <- levels(factor(def.exc))
# reverse <- T
# lst.def <- def.exc == de.lev[1]
# x.def <- x[lst.def]
# y.def <- y[lst.def]
# yy.def <- matrix(, length(x.def), nth)
# pipi.tmp <- numeric(sum(lst.def))
# for (i in (1:nth)) {
# tmp <- pipi[, i]
# tmp <- tmp[lst.def]
# pipi.tmp <- cbind.data.frame(pipi.tmp, tmp)
# }
# pipi.def <- pipi.tmp[, 2:(nth + 1)]
# ycum.def <- matrix(, length(x.def), nth)
# picum.def <- matrix(, length(x.def), nth)
# for (k in 1:nth) {
# yy.def[, k] <- 1 * (y.def == k)
# if (reverse) {
# yy.def[, k] <- rev(yy.def[, k])
# pipi.def[, k] <- rev(pipi.def[, k])
# }
# ycum.def[, k] <- cumsum(yy.def[, k])
# picum.def[, k] <- cumsum(pipi.def[, k])
# if (reverse) {
# ycum.def[, k] <- rev(ycum.def[, k])
# picum.def[, k] <- rev(picum.def[, k])
# }
# }
# lst.exc <- def.exc == de.lev[2]
# x.exc <- x[lst.exc]
# y.exc <- y[lst.exc]
# yy.exc <- matrix(, length(x.exc), nth)
# pipi.tmp <- numeric(sum(lst.exc))
# for (i in (1:nth)) {
# tmp <- pipi[, i]
# tmp <- tmp[lst.exc]
# pipi.tmp <- cbind.data.frame(pipi.tmp, tmp)
# }
# pipi.exc <- pipi.tmp[, 2:(nth + 1)]
# ycum.exc <- matrix(, length(x.exc), nth)
# picum.exc <- matrix(, length(x.exc), nth)
# for (k in 1:nth) {
# yy.exc[, k] <- 1 * (y.exc == k)
# ycum.exc[, k] <- cumsum(yy.exc[, k])
# picum.exc[, k] <- cumsum(pipi.exc[, k])
# }
# x <- c(x.def, x.exc)
# ycum <- rbind(ycum.def, ycum.exc)
# picum <- rbind(picum.def, picum.exc)
# } else {
yy <- matrix(, length(x), nth)
ycum <- matrix(, length(x), nth)
picum <- matrix(, length(x), nth)
for (k in 1:nth) {
yy[, k] <- 1 * (y == k)
ycum[, k] <- cumsum(yy[, k])
picum[, k] <- cumsum(pipi[, k])
}
# }
y.max <- max(ycum)
if (categ.ans < 2) {
if (first.plot == T) {
y.lim <- c(0, y.max)
plot(x, ycum[, 1], xlim = x.lim, ylim = y.lim,
xlab = x.leg, ylab = y.leg, col = 1, cex = cex.tmp,
type = "n")
}
for (k in 1:nth) {
points(x, ycum[, k], pch = paste(k), cex = cex.tmp,
col = col.fct)
lines(x, picum[, k], col = col.fct, lty = k)
}
}
# if (categ.ans == 2) {
# f.graph.window(nth)
# for (ii in 1:nth) {
# plot(x, ycum[, ii], ylim = c(0, max(ycum[,
# ii], picum[, ii])), xlab = x.leg, ylab = y.leg,
# type = "n", col = 1, cex = 0.6)
# points(x, ycum[, ii], pch = "*", cex = 1.5,
# col = 1)
# lines(x, picum[, ii], col = 1, lty = k)
# title(paste("\ncategory ", ii), cex = 0.6)
# }
# first.plot <- F
# }
x.max <- x.high
# if (plot.type == 6) x.max <- 1.1 * (x.high)
if (categ.ans < 2) first.plot <- F
},
# 9
'9' = {
zz <- f.expect.cat(model.type = 1, model.ans, xline,
regr.par, th.par, latent = T, ces.ans = ces.ans,
dtype = dtype, twice = twice, track = track)
y.leg <- paste(y.leg, " (log(latent var))")
upper.uu <- max(upper.uu, exp(max(zz)))
if (sum(is.na(zz)) > 0) return()
if (first.plot == T)
plot(xline, zz, xlim = x.lim,
ylim = c(min(zz), logb(upper.uu)), xlab = x.leg,
ylab = y.leg, type = "n", col = 1)
lines(xline, zz, col = col.fct, lty = 1)
y.max <- max(zz)
x.max <- x.high
first.plot <- F
})
xy.lim.tmp <- c(dum.contr, NA, x.max, 0, y.max, first.plot)
if (track) print("f.lines.tmp in f.lines.cat: END")
return(xy.lim.tmp)
}
nth <- length(th.par)
nr.aa <- max(fct1)
nr.bb <- max(fct2)
regr.par.matr <- f.pars(ans.all)$regr.par.matr
first.plot <- T
kk <- 0
fct2.lev <- levels(factor(fct2))
x.high <- max(x)
upper.uu <- max(regr.par.matr[, 1])
if (plot.type %in% 5:6) {
# Laure added: sort data according to x
xOld <- x
x <- x[order(xOld)]
y <- y[order(xOld)]
y.max <- 0
for (jj in fct2.lev) {
f1 <- fct1[fct2 == as.numeric(jj)]
f1.lev <- as.numeric(levels(factor(f1)))
for (ii in f1.lev) {
y.tmp <- y[fct1 == ii & fct2 == jj]
yy <- matrix(, length(y.tmp), nth)
ycum <- matrix(, length(y.tmp), nth)
for (k in 1:nth) {
yy[, k] <- 1 * (y.tmp == k)
ycum[, k] <- cumsum(yy[, k])
}
y.max <- max(y.max, max(ycum))
y.max <- max(max(ycum))
}
}
}
if (plot.type %in% 5:6) {
kk <- 0
nr.gr <- length(regr.par.matr[, 1])
for (jj in fct2.lev) {
f1 <- fct1[fct2 == as.numeric(jj)]
f1.lev <- levels(factor(f1))
for (ii in f1.lev) {
kk <- kk + 1
if (plot.type != 5 & plot.type != 6)
x.leg <- paste(x.leg, "\ngroup", kk)
xy.lim.tmp <- f.lines.tmp(
x[fct1 == ii & fct2 == jj], y[fct1 == ii & fct2 == jj],
regr.par.matr[kk, ],
th.par, sig.par, nth, dum.contr, plot.type,
col.fct = color[kk], categ.ans = categ.ans,
first.plot = first.plot, model.ans, model.type,
x.leg, y.leg, x.high, upper.uu = upper.uu,
CES.cat = CES.cat, CES = CES, reverse = reverse,
def.exc = def.exc[fct1 == ii & fct2 == jj],
cex.tmp = cex.tmp, show = show, dtype = dtype,
y.max = y.max, x.lim = x.lim, twice = twice)
first.plot <- xy.lim.tmp[6]
}
if ((nr.gr > 1) & (kk < nr.aa * nr.bb)) {
max.xx <- xy.lim[3]
max.yy <- xy.lim[5]
if (plot.type == 6)
max.xx <- max.xx + (max.xx - log10(dum.contr))/200
else max.xx <- max.xx + (max.xx - min(x))/20
}
}
}else for (jj in fct2.lev) {
f1 <- fct1[fct2 == as.numeric(jj)]
f1.lev <- levels(factor(f1))
for (ii in f1.lev) {
kk <- kk + 1
xy.lim.tmp <- f.lines.tmp(
x = x[fct1 == ii & fct2 == jj], y = y[fct1 == ii & fct2 == jj],
regr.par = regr.par.matr[kk, ], th.par, sig.par, nth, dum.contr, plot.type,
col.fct = color[kk], categ.ans = 0, first.plot = first.plot,
model.ans, model.type, x.leg, y.leg, x.high,
upper.uu = upper.uu, CES.cat = CES.cat, CES = CES,
reverse = reverse, def.exc = NA, #def.exc[fct1 == ii & fct2 == jj],
cex.tmp = cex.tmp, show = show,
dtype = dtype, y.max = y.max, x.lim = x.lim,
twice = twice
)
y.max <- xy.lim.tmp[5]
first.plot <- xy.lim.tmp[6]
if (!is.na(CED.matr[1])) {
f.ced.lines(regr.par, nth, dum.contr, CED.matr[kk, ],
x, response = response.matr[kk, ], plot.type,
upper.uu = upper.uu, col.fct = color[kk]
)
}
}
}
if (track) print("f.lines.cat: END")
return(c(0, y.max))
})
}
alfcrisci/bmdModeling documentation built on May 28, 2019, 12:32 a.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.
#' Wrapper function for plotting the results of proast
#' @param ans.all list, with all results that were obtained during the analysis
#' @param track logical, if TRUE (FALSE by default) print the name of the function which is currently being run
#' @param plot.type, integer, defines the type of plot; for cont choose 1 (y vs x),
#' for dtype 2 or 3 choose one of 3 (y.up vs. x), 5 (cumulative y vs. x);
#' default is 1 for cont and 5 else
#' @return list, updated version of ans.all
#' @export
f.plot.all <- function (ans.all, track = FALSE, plot.type = NA) {
if (track) print("f.plot.all")
# MV added default plot type choice
if(is.na(plot.type)){
if(ans.all$cont){
ans.all$plot.type <- 1
} else {
ans.all$plot.type <- 5
}
} else {
ans.all$plot.type <- plot.type
}
if (is.na(ans.all$xy.lim[1]) || ans.all$xy.lim[1] == 0) {
x <- ans.all$x
dum.contr <- min(x[x > 0])/5
ans.all$xy.lim <- c(dum.contr, min(x), max(x), min(ans.all$y), max(ans.all$y))
}
with(ans.all, {
if (cont) {
ans.all$y.lim <- f.plot.con(ans.all, track = track)
if (model.ans != 11){
f.lines.con(ans.all, track = track)
}
if (!is.na(CED[1]) & model.ans != 16) {
f.cedlines.con(ans.all, track = track)
}
} else {
if (dtype %in% c(4, 6, 84) | (dtype == 2 & plot.type < 3)) {
ans.all <- f.plot.frq(ans.all, track = track)
if (model.ans != 14)
ans.all$gr.txt <- f.lines.frq(ans.all, track = track)
f.cedlines.bin(ans.all, track = track)
} else if(dtype %in% 2:3) {
## in f.mm4.cat, code without combi:
ans.all$combi <- F
ans.all$categ.ans <- 0
# f.plot.all(ans.all)
# ans.all.plt <- ans.all #f.model.bb(ans.all = ans.all)
# Error when plotting lines. In Proast61.3, NOTE: Full model may not be applicable for dtype = 2
if (!(model.ans == 14 & dtype == 2))
ans.all$y.lim <- f.lines.cat(ans.all, track = track)
title(ans.all$model.name)#f.mtext(ans.all)
# in f.cat: code with combi.ans:
## ans.all$combi.ans <- 1
## ans.all$categ.ans <- 0
## ans.all$plot.type <- 5
## f.plot.all(ans.all, sep = F)
## f.graph.window(4) equivalent to:
# par(mfcol = c(2, 2), mar = c(3.2, 3.7, 2, 1.5), cex.main = 1.2,
# cex.sub = 0.9, cex.lab = 0.9, cex = 0.9)
# ans.all.plt <- ans.all
# ans.all.plt$modelname <- ""
# ans.all.plt$plot.type <- 3
# ans.all.plt$ced.lines <- TRUE
# ans.all.plt$main <- ""
# ans.all.plt$y.lim <- f.lines.cat(ans.all.plt, track = track)
# title(model.name, col = color[1], cex.main = 0.9, font.main = 1)
## f.mtext(ans.all.plt, track = track)
# ans.all.plt$plot.type <- 9
# ans.all.plt$modelname <- ""
# ans.all.plt$ced.lines <- F
# ans.all.plt$y.lim <- f.lines.cat(ans.all.plt, track = track)
#
# if (!is.na(CED.matr[1])) {
# th.par.cum <- cumsum(th.par)
# CED.vect <- CED.matr[1, ]
# if (ans.scale == 1) {
# CED.vect <- log10(CED.vect)
# x <- log10(x[x > 0]/sf.x)
# }
# for (ii in 1:length(CED.vect))
# lines(c(min(x), CED.vect[ii]), rep(th.par.cum[ii], 2), lty = 2)
# y.leg.2 <- paste(y.leg, " (latent var.uu)")
# if (dtype == 3 & nr.gr == 1) {
# CED.tmp <- round(CED.vect, rep(4, length(CED.vect)))
# CED.tmp <- CED.tmp[CED.tmp < max(x)]
# CED.left <- c(CED.tmp[1]/100, CED.tmp)
# CED.right <- c(CED.tmp, max(x))
# CED.mid <- (CED.left + CED.right)/2
# if (ans.scale == 1)
# CED.mid <- 10^(CED.mid)
# uu <- f.expect.con(model.ans, CED.mid, regr.par,
# fct1 = 1, fct2 = 1, CES = CES, track = track)
# zz <- logb(uu)
# if (ans.scale == 1)
# CED.mid <- log10(CED.mid)
# for (ii in 1:length(scores.orig))
# text(CED.mid[ii], zz[ii], scores.orig[ii])
# }
# }
# th.par.cum <- cumsum(th.par)
# if (mode(CED.matr) == "NULL")
# CED.vect <- 0
# else if (is.na(CED.matr[1]))
# CED.vect <- 0
# else CED.vect <- CED.matr[1, ]
# ans.all.plt$plot.type <- 5
}
}
if (track)
print("f.plot.all : END")
return(ans.all)
})
}
#' Plot the predicted response values (points) and estimated CIs
#' @param ans.all list, with all results that were obtained during the analysis
#' @param track logical, if TRUE (FALSE by default) print the name of the function which is currently being run
#' @return vector, range of the CI bounds on the predicted response values
#' @export
f.plot.con <- function (ans.all, track = FALSE) {
if (track)
print("f.plot.con")
with(ans.all, {
if (is.null(x.leg))
x.leg <- ""
if (is.null(y.leg))
y.leg <- ""
nr.aa <- max(fct1)
nr.bb <- max(fct2)
nr.var <- max(fct3)
twice <- FALSE
if (sum(fct1 != fct2) == 0)
twice <- TRUE
f.means <- function(xx, yy, dtype, sd2.log, plot.type, nn) {
if (track)
print("f.means within f.plot.con")
mean.y <- yy
if (dtype %in% c(1, 5, 15, 101)){
mean.y <- exp(tapply(log(yy), xx, mean))
} else if (dtype == 26){
mean.y <- tapply(sqrt(yy), xx, mean)^2
} else if (dtype == 25){
mean.y <- tapply(yy, xx, mean)
}
if (plot.type %in% c(1, 2, 5, 10, 11))
mean.y <- mean.y
if (plot.type %in% c(3, 4, 6))
mean.y <- log10(mean.y)
if (plot.type %in% c(7, 8))
mean.y <- sqrt(mean.y)
f.var <- function(dtype, xx, yy, sd2.log, nn) {
if (track)
print("f.var within f.plot.con")
if (dtype %in% c(10, 15, 110, 250, 260)) {
SS <- sum(sd2.log * (nn - 1))
var.within <- SS/sum(nn - 1)
df <- sum(nn - 1)
} else {
if (dtype %in% c(1, 5)) {
yy <- log(yy)
} else if (dtype == 26){
yy <- sqrt(yy)
}
mn <- tapply(yy, xx, mean)
nn <- tapply(yy, xx, length)
Vmn <- rep(mn[1], nn[1])
if (dtype %in% c(1, 5)){
for (ii in 2:length(mn)){
Vmn <- c(Vmn, rep(mn[ii], nn[ii]))
}
}
resid <- yy - Vmn
df <- sum(nn - 1)
var.within <- (sum(nn) - 1) * var(resid)/df
}
out.lst <- list(Vsem = sqrt(var.within/nn), df = df)
return(out.lst)
}
var.out <- f.var(dtype, xx, yy, sd2.log, nn)
Vsem <- var.out$Vsem
df <- var.out$df
Vconf <- qt(0.975, df) * Vsem
if (dtype %in% c(1, 5, 10, 15)) {
conf.L <- mean.y/exp(Vconf)
conf.U <- mean.y * exp(Vconf)
} else if (dtype %in% c(25, 250)) {
conf.L <- mean.y - Vconf
conf.U <- mean.y + Vconf
} else {
conf.L <- (sqrt(mean.y) - Vconf)^2
conf.U <- (sqrt(mean.y) + Vconf)^2
}
y.lim.tmp <- range(c(conf.L, conf.U))
if (plot.type %in% c(3, 4, 6)) {
conf.L <- log10(conf.L)
conf.U <- log10(conf.U)
}
if (plot.type %in% c(7, 8)) {
conf.L <- sqrt(conf.L)
conf.U <- sqrt(conf.U)
}
y.lim.tmp <- c(min(conf.L), max(conf.U))
out.lst <- list(conf.L = conf.L, conf.U = conf.U,
y.lim.tmp = y.lim.tmp, mean.y = mean.y)
return(out.lst)
}
x.plt <- x
y.plt <- yy
switch(as.character(plot.type),
'1' = { #1
x.plt <- x
y.plt <- yy
},
'3' = { #3
x.plt <- x
y.plt <- log10(yy)
},
'5' = { #5
y.plt <- yy
x.plt <- sqrt(x)
if (x.leg != "") xleg <- paste("sqrt-", x.leg, sep = "") else xleg <- ""
})
if (dtype %in% c(10, 110, 250, 260)) {
plt.mns <- 3
} else {
plt.mns <- 1
}
cex.1 <- 0.6
cex.2 <- 1.5
mark <- c(1:2, 4:25, 33:500)
# if (nr.aa == 1 && nr.bb == 1 && nr.var == 1) {
#
# if (max(sp.fact) > 1) {
#
# fct2 <- sp.fact
# nr.bb <- max(fct2)
# cat("\n\nno covariates, but subgroups are plotted distinctly\n\n")
#
# }
#
# }
calculationsPlot <- f.means(xx = x.plt, yy = y.plt, dtype = dtype,
sd2.log = sd2.log, plot.type = plot.type, nn = nn)
yRange <- range(c(y.plt, calculationsPlot$y.lim.tmp))
plot(x.plt, y.plt, main = heading, ylim = yRange,
xlab = x.leg, ylab = y.leg, type = "n", col = color[1])
if (nr.aa == 1 & nr.bb == 1 & nr.var > 1) {
nr.bb <- nr.var
fct2 <- fct3
}
if (nr.aa > 1 & nr.bb > 1 && twice) {
nr.aa <- 1
fct1 <- rep(1, length(x))
}
shift.tmp <- 0
shift.abs <- 0
bbb <- (max(x.plt, na.rm = TRUE) - min(x.plt, na.rm = TRUE))/70
zz <- 0
for (jj in 1:nr.bb) {
for (ii in 1:nr.aa) {
x.part <- x.plt[fct1 == ii & fct2 == jj]
y.part <- y.plt[fct1 == ii & fct2 == jj]
y.tmp <- yy[fct1 == ii & fct2 == jj]
if (length(y.part) > 0) {
zz <- zz + 1
if (zz > 500)
cat("\nATTENTION: number of subgroups too large for plotting\n\n")
type.dum <- "p"
lt.y <- 1
if (model.ans == 11) {
lt.y <- 2
type.dum <- "b"
}
x.part <- x.part + shift.tmp
if (plt.mns < 3) {
points(x.part, y.part, col = color[zz], pch = mark[zz],
cex = cex.1, type = "p", lty = lt.y)
}
mean.x <- NA
# Plot mean
# Note: changed compared to original code, because
# Error in xy.coords(x, y) : 'x' and 'y' lengths differ
if (dtype %in% c(1, 5, 15, 25, 26)) {
mean.x <- as.numeric(tapply(x.part, x.part, mean))
mean.y <- f.means(x.part, y.tmp, dtype, plot.type = plot.type)$mean.y
points(mean.x, mean.y, col = color[zz], pch = mark[zz],
cex = 1, type = type.dum, lty = lt.y)
} else {
mean.x <- x.part
mean.y <- y.part
points(mean.x, mean.y, col = color[zz], pch = mark[zz],
cex = 1, type = type.dum, lty = lt.y)
}
if (dtype %in% c(10, 250, 260)) {
sd2.log.part <- sd2.log[fct1 == ii & fct2 == jj]
nn.part <- nn[fct1 == ii & fct2 == jj]
}
out.lst <- f.means(x.part, y.tmp, dtype, sd2.log = sd2.log.part,
plot.type, nn = nn.part)
conf.L.part <- out.lst$conf.L
conf.U.part <- out.lst$conf.U
# Plot confidence limits
for (qq in 1:length(conf.L.part)) {
mean.x[qq] <- mean.x[qq]
lines(rep(mean.x[qq], 2), c(conf.L.part[qq],
conf.U.part[qq]), col = color[zz])
lines(c(mean.x[qq] - bbb, mean.x[qq] +
bbb), rep(conf.L.part[qq], 2), col = color[zz])
lines(c(mean.x[qq] - bbb, mean.x[qq] +
bbb), rep(conf.U.part[qq], 2), col = color[zz])
}
}
shift.tmp <- shift.tmp + shift.abs
}
}
ans.all$y.lim <- range(y.plt)
ans.all$y.lim.CI <- range(y.plt)
if (track)
print("f.plot.con: END")
return(calculationsPlot$y.lim.tmp)
})
}
#' Plot dotted lines for the determined CED value for continous model
#' @param ans.all list, with all results that were obtained during the analysis
#' @param track logical, if TRUE (FALSE by default) print the name of the function which is currently being run
#' @return NULL
#' @export
f.cedlines.con <- function (ans.all, track = FALSE) {
if (track) print("f.cedlines.con")
with(ans.all, {
nr.subgr <- nrow(regr.par.matr)
dum.contr <- xy.lim[1]
# MV changed this
# low.y <- min(y.lim)
low.y <- xy.lim[4]
for (jj in 1:nr.subgr) {
CED.0 <- unlist(CED[jj])
CES.0 <- CES
ES.0 <- f.expect.con(model.ans, CED.0, regr.par.matr[jj,],
fct1 = 1, fct2 = 1, fct5 = 1, CES = CES.0,
increase = increase, track = track)
ES.y <- rep(ES.0, 2)
ES.x <- c(min(x), CED.0)
CED.y <- c(low.y, ES.0)
CED.x <- rep(CED.0, 2)
lines(ES.x, ES.y, lty = 2, lwd = 1)
lines(CED.x, CED.y, lty = 2, lwd = 1)
}
if (track) print("f.cedlines.con: END")
return(NULL)
})
}
#' Plot the curve of the estimated model
#' @param ans.all list, with all results that were obtained during the analysis
#' @param track logical, if TRUE (FALSE by default) print the name of the function which is currently being run
#' @return NULL
#' @export
f.lines.con <- function (ans.all, track = FALSE) {
if (track)
print("f.lines.con")
with(ans.all, {
f.lines.tmp <- function(ans.all.tmp) {
with(ans.all.tmp, {
nbins <- 1000
nbins.0 <- nbins/10
xline <- seq(from = min(x, na.rm = TRUE),
to = max(x, na.rm = TRUE), length = nbins)
twice <- FALSE
expect <- f.expect.con(model.ans, xline, regr.par,
fct1 = rep(1, length(xline)), fct2 = rep(1,
length(xline)), fct3 = 1, CES = CES, twice = twice,
increase = increase, x.mn = x.mn, track = track)
lines(xline, expect, lty = 1, col = col.tmp)
return(NULL)
})
}
ans.all.tmp <- ans.all
ans.all.tmp$dum.contr <- xy.lim[1]
ans.all.tmp$min.x <- xy.lim[2]
ans.all.tmp$max.x <- xy.lim[3]
if ((nr.aa > 1 || nr.bb > 1) && (nr.cc > 1 || nr.dd > 1) && !twice) {
warning("Combination of covariate on c or d and yet another covariate
\ndoes not allow plotting of all individual curves")
}
if (ans.plt == 1)
regr.par.matr <- matrix(regr.par, nrow = 1)
else regr.par.matr <- f.pars(ans.all.tmp)$regr.par.matr
count <- 0
if (nr.aa == 1 && nr.bb == 1 && nr.var > 1){
for (ii in 1:nr.var) {
count <- count + 1
ans.all.tmp$regr.par <- regr.par.matr[1, ]
ans.all.tmp$col.tmp <- color[count]
f.lines.tmp(ans.all.tmp)
}
} else {
kk.max <- nrow(regr.par.matr)
par.bb.old <- NA
par.bb.new <- 0
for (kk in 1:kk.max) {
ans.all.tmp$regr.par <- regr.par.matr[kk, ]
if (model.ans > 1)
par.bb.new <- ans.all.tmp$regr.par[2]
if (ans.plt == 1) {
if (!identical(par.bb.new, par.bb.old) || (model.ans ==
1)) {
count <- count + 1
ans.all.tmp$col.tmp <- color[count]
}
}
else {
count <- count + 1
ans.all.tmp$col.tmp <- color[count]
}
f.lines.tmp(ans.all.tmp)
par.bb.old <- par.bb.new
}
}
if (track)
print("f.lines.con : END")
return(NULL)
})
}
#' Plot the predicted response values (points) and estimated CIs
#' @param ans.all list, with all results that were obtained during the analysis
#' @param track logical, if TRUE (FALSE by default) print the name of the function which is currently being run
#' @return updated ans.all
#' @export
f.plot.frq <- function (ans.all, track = FALSE){
if(track) print("f.plot.frq")
with(ans.all, {
if (length(fct1) == 1)
fct1 <- rep(1, length(x))
if (length(fct2) == 1)
fct2 <- rep(1, length(x))
dum.contr <- xy.lim[1]
y.lim <- c(xy.lim[4], xy.lim[5])
x.lim <- c(xy.lim[2], xy.lim[3])
if (down) y <- 1 - y
if (max(sp.fact) > 1) fct2 <- sp.fact
if (max(fct1) > 1 & max(fct2) > 1 & (sum(fct1 == fct2) == length(fct1)))
twice <- T
shift.abs <- 0
shift.tmp <- 0
mark <- c(1:25, 32:255)
x.lim.plt <- x.lim
y.lim.plt <- y.lim
if (shift > 0) shift.abs <- (max(x) - min(x))/shift
fct1 <- as.numeric(fct1)
fct2 <- as.numeric(fct2)
dum <- max(fct1, fct2)
plot(x, y, #ylim = y.lim, xlim = x.lim, #removed because error if NA
main = heading, xlab = x.leg, ylab = y.leg, type = "n",
col = color[1])
if (shift == 0)
bbb <- (max(x) - min(x))/70 else
bbb <- (max(x) - min(x))/shift
if (max(as.numeric(fct3)) > 1) fct1 <- fct3
zz <- 0
for (jj in 1:max(fct2)) {
for (ii in 1:max(fct1)) {
x.tmp <- x[fct1 == ii & fct2 == jj]
if (!is.na(x.tmp[1])) {
zz <- zz + 1
x.tmp <- x.tmp + shift.tmp
y.tmp <- y[fct1 == ii & fct2 == jj]
n.tmp <- nn[fct1 == ii & fct2 == jj]
if (is.na(cex.1)) {
if (dtype == 6)
cex.1 <- 0.6
else cex.1 <- 1.2
}
points(x.tmp, y.tmp, col = color[zz], pch = mark[zz], cex = cex.1)
if (dtype == 4 & model.ans == 14)
lines(x.tmp, y.tmp, col = color[zz], lty = 2)
if (CI.plt) {
k.tmp <- kk[fct1 == ii & fct2 == jj]
k.tmp <- k.tmp[!is.na(k.tmp)]
n.tmp <- n.tmp[!is.na(n.tmp)]
L025 <- numeric()
L975 <- numeric()
for (qq in 1:length(k.tmp)) {
LL.tmp <- f.LL.bin(k.tmp[qq], n.tmp[qq])
L025[qq] <- LL.tmp[1]
L975[qq] <- LL.tmp[2]
}
if (dtype %in% c(4, 6, 84)) {
points(x.tmp, L025, pch = "-", cex = 1.2, col = color[zz])
points(x.tmp, L975, pch = "-", cex = 1.2, col = color[zz])
for (qq in 1:length(L025)) {
lines(rep(x.tmp[qq], 2), c(L025[qq], L975[qq]), col = color[zz])
lines(c(x.tmp[qq] - bbb, x.tmp[qq] + bbb), rep(L025[qq], 2), col = color[qq])
lines(c(x.tmp[qq] - bbb, x.tmp[qq] + bbb), rep(L975[qq], 2), col = color[zz])
}
}
shift.tmp <- shift.tmp + shift.abs
}
}
}
}
if (!is.na(fct1.full[1]))
if (dtype == 6) {
if (is.null(ans.all$cex.2))
cex.2 <- 1.5
zz <- 0
for (jj in 1:max(fct2.full)) {
for (ii in 1:max(fct1.full)) {
x.tmp <- x.full[fct1.full == ii & fct2.full ==
jj]
if (!is.na(x.tmp[1])) {
zz <- zz + 1
x.tmp <- x.tmp + shift.tmp
y.tmp <- pi.full[fct1.full == ii & fct2.full ==
jj]
if (sum(is.na(pi.full)) == 0) {
if (model.ans == 14)
type.dum <- "b"
else type.dum <- "p"
points(x.tmp, y.tmp, cex = cex.2, pch = mark[zz],
col = color[zz], type = type.dum, lty = 2)
}
}
}
}
}
if (track) print("plot.frq: END")
ans.all$y.lim <- y.lim
ans.all$y.lim.plt <- y.lim.plt
ans.all$shift.abs <- shift.abs
return(ans.all)
})
}
#' Plot the lines
#' @param ans.all list, with all results that were obtained during the analysis
#' @param track logical, if TRUE (FALSE by default) print the name of the function which is currently being run
#' @return ans.all$gr.txt
#' @export
f.lines.frq <- function (ans.all, track = FALSE){
if (track) print("f.lines.frq")
with(ans.all, {
shift.abs <- 0
if (model.ans == 14) return(invisible())
x.high <- max(x)
f.lines.tmp <- function(model.ans, x, regr.par,
th.par, sig.par, plot.type, dtype, l.ty = l.ty,
color, x.high, shift, CES = CES,
twice, pi.full, down) {
if (track) print("f.lines.tmp in f.lines.frq")
sat.mod <- FALSE
if (model.ans == 14) sat.mod <- TRUE
if (!is.na(pi.full[1]) & sat.mod) {
if (dtype %in% c(4, 84))
xline <- x
f(dtype == 6)
xline <- tapply(x, x, mean)
expect <- pi.full
expect <- expect[order(xline)]
xline <- sort(xline)
} else {
xline <- seq(from = min(x), to = x.high, length = 1000)
expect <- f.expect.bin(model.ans, x = xline, regr.par,
CES = CES, ttt = ttt, twice = twice, ces.ans = ces.ans,
CES1 = CES1, CES2 = CES2)
expectdum <- f.expect.bin(model.ans, x = 0, regr.par,
CES = CES, ttt = ttt, twice = twice, ces.ans = ces.ans,
CES1 = CES1, CES2 = CES2)
if (down) {
expect <- 1 - expect
expectdum <- 1 - expectdum
}
pipi <- f.expect.bin(model.ans, x = x, regr.par,
CES = CES, ttt = ttt, twice = twice, ces.ans = ces.ans,
CES1 = CES1, CES2 = CES2)
if (sum(is.na(expect)) > 0) return()
xline <- xline + shift
}
lines(xline, expect, lty = l.ty, col = color)
}
if (max(as.numeric(fct3)) > 1) {
ans.all <- f.pars.frq(ans.all, track = track)
regr.par.matr <- ans.all$regr.par.matr
th.par.vec <- ans.all$th.par.vec
kk.max <- length(th.par.vec)
# MV added this to prevent errors for model.ans = 1 with covariates
if(kk.max != nrow(regr.par.matr))
regr.par.matr <- t(regr.par.matr)
} else {
ans.all <- f.pars(ans.all, track = track)
regr.par.matr <- ans.all$regr.par.matr
# MV changed this to prevent showing too many lines
# kk.max <- length(regr.par.matr[, 1])
kk.max <- nrow(unique(regr.par.matr))
th.par.vec <- rep(th.par, kk.max)
}
shift.tmp <- 0
for (kk in 1:kk.max) {
f.lines.tmp(model.ans, x, regr.par.matr[kk, ],
th.par = th.par.vec[kk], sig.par, plot.type,
dtype = dtype, l.ty = l.ty,
col = color[kk], x.high = x.high,
shift = shift.tmp, CES = CES,
twice = twice, pi.full = pi.full,
down = down)
shift.tmp <- shift.tmp + shift.abs
}
if (track) print("f.lines.frq : END")
return(ans.all$gr.txt)
})
}
#' ?
#' @param ans.all list, with all results that were obtained during the analysis
#' @param track logical, if TRUE (FALSE by default) print the name of the function which is currently being run
#' @return updated ans.all
#' @export
f.pars.frq <- function (ans.all, track = FALSE)
{
if (track) print("f.pars.frq")
with(ans.all, {
if (max(fct1) == 1 & max(fct2) == 1 & max(as.numeric(fct3)) == 1)
regr.par.matr <- matrix(regr.par, nrow = 1)
else {
if (max(as.numeric(fct3)) == 1)
pars <- c(regr.par, th.par)
else {
th.lvm <- f.th.lvm(th.par, fct3, fct3.ref, track = track)
pars <- c(regr.par, th.lvm)
fct3 <- as.numeric(factor(fct3))
}
nrp <- length(regr.par)
regr.par.matr <- numeric()
nr.aa <- max(fct1)
nr.bb <- max(fct2)
nr.th <- max(as.numeric(fct3))
if (!covar.dd) {
if ((model.ans == 11) & (nr.th > 1))
fct1 <- fct3
par.tmp <- rep(NA, length(pars))
if (nrp > (nr.aa + nr.bb))
par.rest <- pars[(nr.aa + nr.bb + 1):nrp]
else par.rest <- numeric()
zz <- 1
for (jj in 1:nr.bb) {
f1 <- fct1[fct2 == jj]
f1.lev <- levels(factor(f1))
for (ii.index in f1.lev) {
ii <- as.numeric(ii.index)
if (max(fct3) > 1) {
f3 <- fct3[fct1 == ii & fct2 == jj]
f3.lev <- levels(factor(f3))
for (kk.index in f3.lev) {
kk <- as.numeric(kk.index)
if (!twice)
gr.txt[zz] <- paste(fct2.txt[jj], fct3.txt[kk], sep = "-")
if (twice)
gr.txt[zz] <- fct2.txt[jj]
par.tmp <- c(pars[ii], pars[nr.aa + jj],
par.rest, pars[nr.aa + nr.bb + length(par.rest) + kk])
regr.par.matr <- rbind(regr.par.matr, par.tmp)
zz <- zz + 1
}
}
if (max(fct3) == 1) {
par.tmp <- c(pars[ii], pars[nr.aa + jj], par.rest, th.par)
regr.par.matr <- rbind(regr.par.matr, par.tmp)
if (!twice)
gr.txt[kk] <- paste(fct1.txt[ii], fct2.txt[jj],
sep = "-")
if (twice)
gr.txt[kk] <- fct1.txt[ii]
}
}
}
n.col <- nrp - nr.aa - nr.bb + 2 + 1
regr.par.matr <- matrix(regr.par.matr[, 1:n.col],
ncol = n.col)
}
if (covar.dd) {
cc <- pars[nr.aa + nr.bb + 1]
for (jj in 1:nr.bb) {
par.tmp <- c(pars[jj], pars[nr.aa + jj], cc,
pars[nr.aa + nr.bb + 1 + jj])
regr.par.matr <- rbind(regr.par.matr, par.tmp)
}
}
}
npar <- length(regr.par.matr[1, ])
th.par.vec <- regr.par.matr[, npar]
regr.par.matr <- regr.par.matr[, 1:(npar - 1)]
# MV added this to prevent errors
if(!is.matrix(regr.par.matr))
regr.par.matr <- matrix(regr.par.matr, nrow = 1)
ans.all$regr.par.matr <- regr.par.matr
ans.all$gr.txt <- gr.txt
ans.all$th.par.vec <- th.par.vec
if (track) print("f.pars.frq: END")
return(ans.all)
})
}
#' compute th.lvm?
#' @param th.par parameter theta
#' @param fct3 initial value for theta
#' @param fct3.ref theta reference?
#' @param track logical, if TRUE (FALSE by default) print the name of the function which is currently being run
#' @return th.lvm
#' @export
f.th.lvm <- function (th.par, fct3, fct3.ref, track = FALSE) {
if (track) print("f.th.lvm")
nr.th <- length(levels(as.factor(fct3)))
fct3.lev <- as.numeric(levels(as.factor(fct3)))
th.lvm <- rep(th.par, nr.th)
th.lvm <- th.lvm - (log(fct3.ref) - log(fct3.lev))
if (track) print("f.th.lvm: END")
return(th.lvm)
}
#' Plot dotted lines for the determined CED value for categorical model
#' @param ans.all list, with all results that were obtained during the analysis
#' @param inputResponse numeric, value for the response at estimated CED;
#' default value is NA
#' @param inputCED numeric, value for the CED; default value is NA
#' @param track logical, if TRUE (FALSE by default) print the name of the function which is currently being run
#' @return NULL
#' @export
f.cedlines.bin <- function (ans.all, inputResponse = NA, inputCED = NA,
track = FALSE) {
if (track) print("f.cedlines.bin")
with(ans.all, {
# MV added these conditions to use plot for model averaging
if(is.na(inputResponse))
response <- response.matr else
response <- inputResponse
if(is.na(inputCED))
CED <- CED.matr else
CED <- inputCED
dum.contr <- xy.lim[1]
nr.lev <- max(max(fct1, fct2))
ES.y <- matrix(response[1], ncol = 2)
ES.x <- matrix(c(x[1], CED[1]), ncol = 2)
if (nr.lev > 1)
for (ii in (2:nr.lev)) {
ES.y <- rbind(ES.y, rep(response[ii], 2))
ES.x <- rbind(ES.x, c(x[1], CED[ii]))
}
CED.y <- matrix(c(min(y), response[1]), ncol = 2)
CED.x <- matrix(rep(CED[1], 2), ncol = 2)
if (nr.lev > 1)
for (ii in (2:nr.lev)) {
CED.y <- rbind(CED.y, c(min(y), response[ii]))
CED.x <- rbind(CED.x, rep(CED[ii], 2))
}
for (ii in (1:nr.lev)) {
if (down) {
ES.y[ii, ] <- 1 - ES.y[ii, ]
CED.y[ii, 2] <- 1 - CED.y[ii, 2]
}
lines(ES.x[ii, ], ES.y[ii, ], lty = 2)
lines(CED.x[ii, ], CED.y[ii, ], lty = 2)
}
return(NULL)
})
}
#' Calculate confidence intervals for each of the plotted points for quantal response
#' @param k integer, number of events in given group
#' @param n integer, number of observations in given group
#' @param n.pi ?
#' @param conf ?
#' @param step ?
#' @return vector with lower and upper limit
#' @export
f.LL.bin <- function (k = 4, n = 20, n.pi = 100, conf = 0.9, step = 1.2)
{
compl <- 0
alfa <- (1 - conf)/2
if ((k/n) > 0.5) {
compl <- T
k <- n - k
}
pi.min <- (k/n)
pi.max <- 1
pi.range <- 0.999 * seq(pi.min, pi.max, (pi.max - pi.min)/n.pi)
prob1 <- pbinom(k, n, pi.range)
if (sum(pi.range < 0 | pi.range > 1) > 0) {
}
dum1 <- pi.range[prob1 <= alfa]
while (length(dum1) > n.pi/2) {
pi.min <- pi.min/step
pi.max <- pi.max/step
pi.range <- 0.999 * seq(pi.min, pi.max, (pi.max - pi.min)/n.pi)
prob1 <- pbinom(k, n, pi.range)
if (sum(pi.range < 0 | pi.range > 1) > 0) {
}
dum1 <- pi.range[prob1 <= alfa]
}
dum2 <- pi.range[prob1 >= alfa]
P2 <- (max(prob1[prob1 <= alfa]))
P1 <- (min(prob1[prob1 >= alfa]))
L.up <- ((P1 - alfa) * (min(dum1) - max(dum2)))/(P1 - P2) +
max(dum2)
pi.min <- (k/n)
pi.max <- 1
pi.range <- 0.999 * seq(pi.min, pi.max, (pi.max - pi.min)/n.pi)
if (k == 0)
L.low <- 0
else {
prob2 <- 1 - pbinom(k - 1, n, pi.range)
if (sum(pi.range < 0 | pi.range > 1) > 0) {
}
if (sum(pi.range < 0 | pi.range > 1) > 0) {
}
dum1 <- pi.range[prob2 <= alfa]
while (length(dum1) < n.pi/2) {
pi.max <- pi.max/step
pi.min <- pi.min/step
pi.range <- 0.999 * seq(pi.min, pi.max, (pi.max -
pi.min)/n.pi)
prob2 <- 1 - pbinom((k - 1), n, pi.range)
if (sum(pi.range < 0 | pi.range > 1) > 0) {
}
dum1 <- pi.range[prob2 <= alfa]
}
dum2 <- pi.range[prob2 >= alfa]
P2 <- (max(prob2[prob2 <= alfa]))
P1 <- (min(prob2[prob2 >= alfa]))
L.low <- min(dum2) - ((P1 - alfa) * (min(dum2) - max(dum1)))/(P1 -
P2)
}
if (compl) {
L.low.tmp <- L.low
L.low <- 1 - L.up
L.up <- 1 - L.low.tmp
}
return(c(L.low, L.up))
}
#' plot lines for categorical model
#'
#' Note: used in the f.plot.all function, for dtype == 2, 3,
#' and where plot.type is set to 3, 5, 9?
#' so only this part of the code is retained
#' @param ans.all list, with all results that were obtained during the analysis
#' @param track logical, if TRUE (FALSE by default) print the name of the function which is currently being run
#' @return c(0, y.max)
#' @export
f.lines.cat <- function (ans.all, track = FALSE){
if (track) print("f.lines.cat")
if (length(ans.all$cex.tmp) == 0) cex.tmp <- 1
with(ans.all, {
dum.contr <- xy.lim[1]
x.lim <- c(xy.lim[2], xy.lim[3])
y.max <- xy.lim[5]
f.ced.lines <- function(regr.par, nth, dum.contr, CED,
x, response, plot.type, upper.uu = upper.uu, col.fct) {
if (track) print("f.ced.lines in f.lines.cat")
switch(as.character(plot.type),
'3' = for (k in 1:nth)
lines(rep(CED[k], 2), c(0, response[k]),
lty = 2, col = col.fct),
'9' = for (k in 1:nth)
lines(rep(CED[k], 2), c(min(logb(regr.par[2])), th.par[k]),
lty = 2, col = col.fct)
)
switch(as.character(plot.type),
# 3
'3' = for (k in 1:nth)
lines(c(0, CED[k]), c(response[k], response[k]),
lty = 2, col = col.fct),
# 9
'9' = for (k in 1:nth)
lines(c(0, CED[k]), rep(th.par[k], 2), lty = 2, col = col.fct)
)
# nothing is done for plot.type == 5
}
f.lines.tmp <- function(x, y, regr.par, th.par, sig.par,
nth, dum.contr, plot.type, col.fct, categ.ans = categ.ans,
first.plot, model.ans, model.type, x.leg, y.leg,
x.high, upper.uu = upper.uu, CES.cat, CES, reverse = reverse,
def.exc, cex.tmp, show, dtype, y.max, x.lim, twice) {
cat.txt <- c("category 1", "category 2", "category 3",
"category 4", "category 5", "category 6", "category 7",
"category 8")
if (track) print("f.lines.tmp in f.lines.cat")
xline <- seq(from = dum.contr, to = x.high, length = 100)
if (model.type == 1) {
#MV error when plot for dtype = 2, so added first part using dtype = 2 instead of default dtype = 6
if(dtype == 2){
pipi <- matrix(f.expect.bin(model.ans, x, regr.par,
CES = CES, ces.ans = ces.ans, dtype = dtype, track = track), ncol = 1)
piline <- matrix(f.expect.bin(model.ans, xline,
regr.par, CES = CES, ces.ans = ces.ans, dtype = dtype, track = track), ncol = 1)
} else {
pipi <- matrix(f.expect.bin(model.ans, x, regr.par,
CES = CES, ces.ans = ces.ans, track = track), ncol = 1)
piline <- matrix(f.expect.bin(model.ans, xline,
regr.par, CES = CES, ces.ans = ces.ans, track = track), ncol = 1)
}
} else {
pipi <- f.expect.cat(model.type, model.ans, x,
regr.par, th.par, sig.par, CES = CES, CES.cat = CES.cat,
ces.ans = ces.ans, dtype = dtype, twice = twice, track = track)
piline <- f.expect.cat(model.type, model.ans,
xline, regr.par, th.par, sig.par, CES = CES,
CES.cat = CES.cat, ces.ans = ces.ans, dtype = dtype,
twice = twice, track = track)
}
switch(as.character(plot.type),
# 3
'3' = {
pi.up <- matrix(0, length(xline), (nth + 1))
if (first.plot == T)
plot(xline, pi.up[, 1],
xlim = x.lim, ylim = c(0, 1), xlab = x.leg,
ylab = paste(y.leg, "(fraction of scores)"),
type = "n", col = col.fct)
for (k in 1:nth) {
dum <- nth + 1 - k
pi.up[, dum] <- piline[, dum] + pi.up[, (dum + 1)]
}
for (k in 1:nth)
lines(xline, pi.up[, k], col = col.fct, lty = k)
y.max <- 1
x.max <- x.high
first.plot <- F
},
'5' = {
y.leg <- paste(y.leg, " (cumu)")
# if (!is.na(def.exc[1])) {
# de.lev <- levels(factor(def.exc))
# reverse <- T
# lst.def <- def.exc == de.lev[1]
# x.def <- x[lst.def]
# y.def <- y[lst.def]
# yy.def <- matrix(, length(x.def), nth)
# pipi.tmp <- numeric(sum(lst.def))
# for (i in (1:nth)) {
# tmp <- pipi[, i]
# tmp <- tmp[lst.def]
# pipi.tmp <- cbind.data.frame(pipi.tmp, tmp)
# }
# pipi.def <- pipi.tmp[, 2:(nth + 1)]
# ycum.def <- matrix(, length(x.def), nth)
# picum.def <- matrix(, length(x.def), nth)
# for (k in 1:nth) {
# yy.def[, k] <- 1 * (y.def == k)
# if (reverse) {
# yy.def[, k] <- rev(yy.def[, k])
# pipi.def[, k] <- rev(pipi.def[, k])
# }
# ycum.def[, k] <- cumsum(yy.def[, k])
# picum.def[, k] <- cumsum(pipi.def[, k])
# if (reverse) {
# ycum.def[, k] <- rev(ycum.def[, k])
# picum.def[, k] <- rev(picum.def[, k])
# }
# }
# lst.exc <- def.exc == de.lev[2]
# x.exc <- x[lst.exc]
# y.exc <- y[lst.exc]
# yy.exc <- matrix(, length(x.exc), nth)
# pipi.tmp <- numeric(sum(lst.exc))
# for (i in (1:nth)) {
# tmp <- pipi[, i]
# tmp <- tmp[lst.exc]
# pipi.tmp <- cbind.data.frame(pipi.tmp, tmp)
# }
# pipi.exc <- pipi.tmp[, 2:(nth + 1)]
# ycum.exc <- matrix(, length(x.exc), nth)
# picum.exc <- matrix(, length(x.exc), nth)
# for (k in 1:nth) {
# yy.exc[, k] <- 1 * (y.exc == k)
# ycum.exc[, k] <- cumsum(yy.exc[, k])
# picum.exc[, k] <- cumsum(pipi.exc[, k])
# }
# x <- c(x.def, x.exc)
# ycum <- rbind(ycum.def, ycum.exc)
# picum <- rbind(picum.def, picum.exc)
# } else {
yy <- matrix(, length(x), nth)
ycum <- matrix(, length(x), nth)
picum <- matrix(, length(x), nth)
for (k in 1:nth) {
yy[, k] <- 1 * (y == k)
ycum[, k] <- cumsum(yy[, k])
picum[, k] <- cumsum(pipi[, k])
}
# }
y.max <- max(ycum)
if (categ.ans < 2) {
if (first.plot == T) {
y.lim <- c(0, y.max)
plot(x, ycum[, 1], xlim = x.lim, ylim = y.lim,
xlab = x.leg, ylab = y.leg, col = 1, cex = cex.tmp,
type = "n")
}
for (k in 1:nth) {
points(x, ycum[, k], pch = paste(k), cex = cex.tmp,
col = col.fct)
lines(x, picum[, k], col = col.fct, lty = k)
}
}
# if (categ.ans == 2) {
# f.graph.window(nth)
# for (ii in 1:nth) {
# plot(x, ycum[, ii], ylim = c(0, max(ycum[,
# ii], picum[, ii])), xlab = x.leg, ylab = y.leg,
# type = "n", col = 1, cex = 0.6)
# points(x, ycum[, ii], pch = "*", cex = 1.5,
# col = 1)
# lines(x, picum[, ii], col = 1, lty = k)
# title(paste("\ncategory ", ii), cex = 0.6)
# }
# first.plot <- F
# }
x.max <- x.high
# if (plot.type == 6) x.max <- 1.1 * (x.high)
if (categ.ans < 2) first.plot <- F
},
# 9
'9' = {
zz <- f.expect.cat(model.type = 1, model.ans, xline,
regr.par, th.par, latent = T, ces.ans = ces.ans,
dtype = dtype, twice = twice, track = track)
y.leg <- paste(y.leg, " (log(latent var))")
upper.uu <- max(upper.uu, exp(max(zz)))
if (sum(is.na(zz)) > 0) return()
if (first.plot == T)
plot(xline, zz, xlim = x.lim,
ylim = c(min(zz), logb(upper.uu)), xlab = x.leg,
ylab = y.leg, type = "n", col = 1)
lines(xline, zz, col = col.fct, lty = 1)
y.max <- max(zz)
x.max <- x.high
first.plot <- F
})
xy.lim.tmp <- c(dum.contr, NA, x.max, 0, y.max, first.plot)
if (track) print("f.lines.tmp in f.lines.cat: END")
return(xy.lim.tmp)
}
nth <- length(th.par)
nr.aa <- max(fct1)
nr.bb <- max(fct2)
regr.par.matr <- f.pars(ans.all)$regr.par.matr
first.plot <- T
kk <- 0
fct2.lev <- levels(factor(fct2))
x.high <- max(x)
upper.uu <- max(regr.par.matr[, 1])
if (plot.type %in% 5:6) {
# Laure added: sort data according to x
xOld <- x
x <- x[order(xOld)]
y <- y[order(xOld)]
y.max <- 0
for (jj in fct2.lev) {
f1 <- fct1[fct2 == as.numeric(jj)]
f1.lev <- as.numeric(levels(factor(f1)))
for (ii in f1.lev) {
y.tmp <- y[fct1 == ii & fct2 == jj]
yy <- matrix(, length(y.tmp), nth)
ycum <- matrix(, length(y.tmp), nth)
for (k in 1:nth) {
yy[, k] <- 1 * (y.tmp == k)
ycum[, k] <- cumsum(yy[, k])
}
y.max <- max(y.max, max(ycum))
y.max <- max(max(ycum))
}
}
}
if (plot.type %in% 5:6) {
kk <- 0
nr.gr <- length(regr.par.matr[, 1])
for (jj in fct2.lev) {
f1 <- fct1[fct2 == as.numeric(jj)]
f1.lev <- levels(factor(f1))
for (ii in f1.lev) {
kk <- kk + 1
if (plot.type != 5 & plot.type != 6)
x.leg <- paste(x.leg, "\ngroup", kk)
xy.lim.tmp <- f.lines.tmp(
x[fct1 == ii & fct2 == jj], y[fct1 == ii & fct2 == jj],
regr.par.matr[kk, ],
th.par, sig.par, nth, dum.contr, plot.type,
col.fct = color[kk], categ.ans = categ.ans,
first.plot = first.plot, model.ans, model.type,
x.leg, y.leg, x.high, upper.uu = upper.uu,
CES.cat = CES.cat, CES = CES, reverse = reverse,
def.exc = def.exc[fct1 == ii & fct2 == jj],
cex.tmp = cex.tmp, show = show, dtype = dtype,
y.max = y.max, x.lim = x.lim, twice = twice)
first.plot <- xy.lim.tmp[6]
}
if ((nr.gr > 1) & (kk < nr.aa * nr.bb)) {
max.xx <- xy.lim[3]
max.yy <- xy.lim[5]
if (plot.type == 6)
max.xx <- max.xx + (max.xx - log10(dum.contr))/200
else max.xx <- max.xx + (max.xx - min(x))/20
}
}
}else for (jj in fct2.lev) {
f1 <- fct1[fct2 == as.numeric(jj)]
f1.lev <- levels(factor(f1))
for (ii in f1.lev) {
kk <- kk + 1
xy.lim.tmp <- f.lines.tmp(
x = x[fct1 == ii & fct2 == jj], y = y[fct1 == ii & fct2 == jj],
regr.par = regr.par.matr[kk, ], th.par, sig.par, nth, dum.contr, plot.type,
col.fct = color[kk], categ.ans = 0, first.plot = first.plot,
model.ans, model.type, x.leg, y.leg, x.high,
upper.uu = upper.uu, CES.cat = CES.cat, CES = CES,
reverse = reverse, def.exc = NA, #def.exc[fct1 == ii & fct2 == jj],
cex.tmp = cex.tmp, show = show,
dtype = dtype, y.max = y.max, x.lim = x.lim,
twice = twice
)
y.max <- xy.lim.tmp[5]
first.plot <- xy.lim.tmp[6]
if (!is.na(CED.matr[1])) {
f.ced.lines(regr.par, nth, dum.contr, CED.matr[kk, ],
x, response = response.matr[kk, ], plot.type,
upper.uu = upper.uu, col.fct = color[kk]
)
}
}
}
if (track) print("f.lines.cat: END")
return(c(0, y.max))
})
}
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.