Last updated on 2020-02-19 10:48:47 CET.
| Flavor | Version | Tinstall | Tcheck | Ttotal | Status | Flags |
|---|---|---|---|---|---|---|
| r-devel-linux-x86_64-debian-clang | 1.0.3 | 16.84 | 122.56 | 139.40 | ERROR | |
| r-devel-linux-x86_64-debian-gcc | 1.0.3 | 15.70 | 94.01 | 109.71 | ERROR | |
| r-devel-linux-x86_64-fedora-clang | 1.0.3 | 174.19 | ERROR | |||
| r-devel-linux-x86_64-fedora-gcc | 1.0.3 | 164.73 | ERROR | |||
| r-devel-windows-ix86+x86_64 | 1.0.3 | 36.00 | 138.00 | 174.00 | OK | |
| r-devel-windows-ix86+x86_64-gcc8 | 1.0.3 | 47.00 | 185.00 | 232.00 | OK | |
| r-patched-linux-x86_64 | 1.0.3 | 13.74 | 114.28 | 128.02 | OK | |
| r-patched-solaris-x86 | 1.0.3 | 238.80 | OK | |||
| r-release-linux-x86_64 | 1.0.3 | 13.04 | 114.15 | 127.19 | OK | |
| r-release-windows-ix86+x86_64 | 1.0.3 | 28.00 | 119.00 | 147.00 | OK | |
| r-release-osx-x86_64 | 1.0.3 | OK | ||||
| r-oldrel-windows-ix86+x86_64 | 1.0.3 | 19.00 | 128.00 | 147.00 | OK | |
| r-oldrel-osx-x86_64 | 1.0.3 | OK |
Version: 1.0.3
Check: examples
Result: ERROR
Running examples in 'bhrcr-Ex.R' failed
The error most likely occurred in:
> base::assign(".ptime", proc.time(), pos = "CheckExEnv")
> ### Name: calculatePCE
> ### Title: WWARN Parasite Clearance Estimator (PCE)
> ### Aliases: calculatePCE
>
> ### ** Examples
>
> ## Don't show:
> data("pursat")
> data = pursat[pursat["id"] <= 80 & pursat["id"] > 70,]
> output <- calculatePCE(data = data, detect.limit = 15, outlier.detect = TRUE)
Calculating WWARN PCE Estimates...
[1] "Replacing a 0 parasite level (at data point 18) with detection level: data set 71"
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
bhrcr
--- call from context ---
lagReg_tobit(R, data_in, data_out, i, start_ind, end_ind, code,
Condition_on_24_hours, fact, Threshold2, Threshold3, Detec_limit,
DT, MaxIters, plotting, Name, ind_DL, FirstHours, MaxDiffInFirstHours,
NeededInFirstHours, TimeRes1, TimeRes2, Threshold4, imageType)
--- call from argument ---
if (class(x_CI) == "try-error") {
x_CI = array(NA, c(2, 2))
} else if (is.na(x_CI[2, 1])) {
x_CI[2, 1] = NA
}
--- R stacktrace ---
where 1: lagReg_tobit(R, data_in, data_out, i, start_ind, end_ind, code,
Condition_on_24_hours, fact, Threshold2, Threshold3, Detec_limit,
DT, MaxIters, plotting, Name, ind_DL, FirstHours, MaxDiffInFirstHours,
NeededInFirstHours, TimeRes1, TimeRes2, Threshold4, imageType)
where 2: calculateParasiteClearance(Name = "pceestimates", minpara = detect.limit,
data1 = data)
where 3: withCallingHandlers(expr, message = function(c) invokeRestart("muffleMessage"))
where 4: suppressMessages(calculateParasiteClearance(Name = "pceestimates",
minpara = detect.limit, data1 = data))
where 5: withCallingHandlers(expr, warning = function(w) invokeRestart("muffleWarning"))
where 6: suppressWarnings(suppressMessages(calculateParasiteClearance(Name = "pceestimates",
minpara = detect.limit, data1 = data)))
where 7: calculatePCE(data = data, detect.limit = 15, outlier.detect = TRUE)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (R, data_in, data_out, i, start_ind, end_ind, code,
Condition_on_24_hours, fact, Threshold2, Threshold3, Detect_limit,
DT, MaxIters, plotting, Name, ind_DL, FirstHours, MaxDiffInFirstHours,
NeededInFirstHours, TimeRes1, TimeRes2, Threshold4, imageType = c("pdf",
"png"))
{
par(ask = FALSE)
imageType <- match.arg(imageType)
indexes = start_ind:end_ind
temp2 = which(R[start_ind:end_ind, 4] == 0)
temp3 = indexes[temp2]
data_out_store = data_out
data_in_store = data_in
data_in = data_in[temp2]
data_out = data_out[temp2]
ldata_out = log(data_out)
n2 = length(ldata_out)
R[temp3, 5] = ldata_out
lag_phase_attempted = 2
MAXREG = 0
filename2 <- sprintf("%s_%s.%s", Name, code, imageType)
NUM = data_in[which(data_in <= FirstHours)]
if (n2 == 3) {
NE = 0
}
else if (length(NUM) <= 1) {
NE = 1
}
else if ((max(diff(NUM)) > MaxDiffInFirstHours) | (length(NUM) <
NeededInFirstHours)) {
NE = 1
}
else {
NE = 0
}
LP_temp = which(data_out != 0)
LLP = LP_temp[length(LP_temp)]
LastPos = data_out[LLP]
k = NULL
tlag = data_in[1]
if (n2 > 3) {
Y = log(data_out[1:(n2 - 1)])
X = data_in[1:(n2 - 1)]
cf = coef(lm(Y ~ X))
x_CI = try(confint(nls(Y ~ beta * (X - x0) - log(DT),
start = c(beta = cf[[2]], x0 = -cf[[1]]/cf[[2]]))))
if (class(x_CI) == "try-error") {
x_CI = array(NA, c(2, 2))
}
else if (is.na(x_CI[2, 1])) {
x_CI[2, 1] = NA
}
}
prof_type = -9999
if (n2 < 3) {
xx <- sprintf("TOBIT: There are not enough time points - no model fitted: data set %s",
code)
R[start_ind, 56] = 1
}
else if (data_out[1] < Threshold2) {
xx <- sprintf("TOBIT: Baseline level is too low - no model fitted: data set %s",
code)
R[start_ind, 56] = 2
}
else if (n2 == 3) {
ldata_out = log(data_out)
R[temp3, 5] = ldata_out
Linear = fitLinearAndStorage(data_in, ldata_out, start_ind,
temp3, R)
glm.linear = Linear$glm.linear
R = Linear$R
AIC1 = Linear$AIC1
mod = glm.linear
LinPart = seq(1, n2)
k = R[start_ind, 21]
alin = mod$coeff[[1]]
R[start_ind, 53] = R[start_ind, 22]
xx <- sprintf("TOBIT: Fitting a linear model since it has the lowest AIC n = 3: data set %s",
code)
prof_type = 1
R[start_ind, 44] = R[start_ind, 23]
}
else if ((data_in[n2] > x_CI[[2, 1]]) & (ldata_out[n2 - 1] >
Threshold3)) {
xx <- sprintf("TOBIT: zero is not informative- no model is fitted: data set %s",
code)
R[start_ind, 56] = 3.1
}
else if (NE == 1) {
ldata_out = log(data_out)
R[temp3, 5] = ldata_out
Linear = try(fitLinearAndStorage_tobit(data_in, ldata_out,
start_ind, temp3, R, DT, MaxIters))
glm.linear = Linear$glm.linear
R = Linear$R
AIC1 = Linear$AIC1
alin = Linear$a1
mod = glm.linear
LinPart = seq(1, n2)
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
k = R[start_ind, 21]
R[start_ind, 53] = R[start_ind, 22]
xx <- sprintf("TOBIT: Fitting a linear model since it has the lowest AIC n = 4: data set %s",
code)
prof_type = 1
R[start_ind, 44] = R[start_ind, 23]
}
else if (n2 == 4) {
ldata_out = log(data_out)
R[temp3, 5] = ldata_out
Linear = try(fitLinearAndStorage_tobit(data_in, ldata_out,
start_ind, temp3, R, DT, MaxIters))
glm.linear = Linear$glm.linear
R = Linear$R
AIC1 = Linear$AIC1
alin = Linear$a1
mod = glm.linear
LinPart = seq(1, n2)
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
k = R[start_ind, 21]
R[start_ind, 53] = R[start_ind, 22]
xx <- sprintf("TOBIT: Fitting a linear model since it has the lowest AIC n = 4: data set %s",
code)
prof_type = 1
R[start_ind, 44] = R[start_ind, 23]
}
else if (n2 == 5) {
Linear = try(fitLinearAndStorage_tobit(data_in, ldata_out,
start_ind, temp3, R, DT, MaxIters))
glm.linear = Linear$glm.linear
R = Linear$R
AIC1 = Linear$AIC1
alin = Linear$a1
Quad = try(fitQuadraticAndStorage_tobit(data_in, ldata_out,
start_ind, temp3, R, DT, MaxIters))
glm.quad = Quad$glm.quad
R = Quad$R
AIC2 = Quad$AIC2
scaledQ = Quad$scaled
ind = which(c(AIC1, AIC2) == min(c(AIC1, AIC2)))
ind = which(c(AIC1, AIC2) == min(c(AIC1, AIC2)))
if (data_out[2] > (1 + fact) * data_out[1]) {
data_in_mod = data_in[2:n2]
ldata_out_mod = ldata_out[2:n2]
tlag = data_in_mod[1]
m1 <- lm(ldata_out_mod ~ poly(data_in_mod, 1, raw = TRUE))
MaxReg <- try(tobit(ldata_out_mod ~ poly(data_in_mod,
1, raw = TRUE), left = log(DT), iter.max = MaxIters))
if (class(MaxReg)[1] == "try-error") {
yy <- sprintf("Fixing error in tobit fitting, data set: %s",
code)
print(yy)
options(warn = 1)
MaxReg <- tobit(ldata_out_mod ~ poly(data_in_mod,
1, raw = TRUE), left = log(DT), init = coef(m1),
iter.max = MaxIters)
options(warn = 2)
}
AICMaxReg = AIC(MaxReg)
ss_1 = sum(residuals(glm.linear)[2:(n2 - 1)]^2)
ss_2 = sum(residuals(glm.quad)[2:(n2 - 1)]^2)
ss_3 = sum(residuals(MaxReg)[1:3]^2)
ind = which(c(ss_1, ss_2, ss_3) == min(c(ss_1, ss_2,
ss_3)))
}
if (ind == 1) {
mod = glm.linear
LinPart = seq(1, n2)
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
tlag = data_in[1]
k = R[start_ind, 21]
R[start_ind, 53] = R[start_ind, 22]
xx <- sprintf("TOBIT: Fitting a linear model since it has the lowest SS n = 5: data set %s",
code)
prof_type = 1
R[start_ind, 44] = R[start_ind, 23]
}
else if (ind == 2) {
mod = glm.quad
mu = mean(data_in)
sd = sd(data_in)
a1 = coef(mod)[[1]]
b1 = coef(mod)[[2]]
c1 = coef(mod)[[3]]
if (scaledQ == 1) {
c = a1 - b1 * mu/sd + c1 * mu^2/sd^2
b = b1/sd - 2 * mu * c1/sd^2
a = c1/sd^2
}
else {
c = coef(mod)[[1]]
b = coef(mod)[[2]]
a = coef(mod)[[3]]
}
mod.fitted <- fitted.values(mod)
slopes = diff(mod.fitted)/diff(data_in)
c_points = c(data_in[1], data_in[length(data_in) -
1])
c_devs = 2 * a * c_points + b
ind_max = 2
fastest_rate = c_devs[ind_max]
ratios = fastest_rate/slopes
LinPart = c(which((ratios <= 5) & (ratios >= 0)),
n2)
ldattemp = mod.fitted[LinPart]
if (a < 0) {
xx <- sprintf("TOBIT: Fitting a quadratic (convex) model since it has the lowest SS and n = 5: data set %s",
code)
prof_type = 3
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
ldattemp = mod.fitted[LinPart]
tlag = data_in[LinPart[1]]
dattemp = data_in[LinPart]
mod <- glm(ldattemp ~ dattemp)
k = mod$coeff[[2]]
alin = mod$coeff[[1]]
}
else if (a >= 0) {
xx <- sprintf("TOBIT: Fitting a linear model since we have identified a concave quadratic (quadratic has the lowest SS) and n = 5: data set %s",
code)
prof_type = 2
LinPart = seq(1, n2)
ldattemp = ldata_out[LinPart]
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
tlag = data_in[LinPart[1]]
mod <- glm.linear
k = R[start_ind, 21]
}
if (tlag == data_in[1]) {
LinPart = seq(1, n2)
LinPart = LinPart[1:(length(LinPart) - 1)]
mod <- glm.linear
k = R[start_ind, 21]
to_include_R2 = LinPart
}
else {
to_include_R2 = seq(1, length(LinPart))
}
temp3 = indexes[temp2]
R[start_ind, 44] = 1 - sum((fitted(mod)[to_include_R2] -
ldata_out[LinPart])^2)/sum((ldata_out[LinPart] -
mean(ldata_out[LinPart]))^2)
SE_temp = R[start_ind, 22]
if (tlag > data_in[1]) {
SE_temp = sqrt(diag(vcov(mod)))[[2]]
}
R[start_ind, 53] = SE_temp
R[start_ind, 50] <- mod.fitted[which(data_in == tlag)[1]] -
mod.fitted[1]
R[start_ind, 50] <- mod.fitted[which(data_in == tlag)[1]] -
mod.fitted[1]
glm.quad.fitted <- fitted.values(glm.quad)
R[start_ind, 51] <- max(glm.quad.fitted) - glm.quad.fitted[which(data_in ==
tlag)[1]]
}
else if (ind == 3) {
mod <- MaxReg
AIC1 = AIC(mod)
LinPart = 1:(n2 - 1)
if (LinPart[length(LinPart)] == (n2 - 1)) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
k = coef(mod)[[2]]
alin = mod$coeff[[1]]
SEs2 = sqrt(diag(vcov(mod)))
R[start_ind, 53] = SEs2[[2]]
xx <- sprintf("TOBIT: Fitting a linear model from max value since we only have 5 data points and an initial increase in parasite level (also - lowest SS): data set %s",
code)
prof_type = 5
R[start_ind, 44] = 1 - sum(residuals(mod)^2)/sum((ldata_out_mod -
mean(ldata_out_mod))^2)
MAXREG = 1
}
}
else {
Linear = fitLinearAndStorage_tobit(data_in, ldata_out,
start_ind, temp3, R, DT, MaxIters)
glm.linear = Linear$glm.linear
R = Linear$R
AIC1 = Linear$AIC1
alin = Linear$a1
Quad = fitQuadraticAndStorage_tobit(data_in, ldata_out,
start_ind, temp3, R, DT, MaxIters)
glm.quad = Quad$glm.quad
R = Quad$R
AIC2 = Quad$AIC2
scaledQ = Quad$scaled
Cubic = fitCubicAndStorage_tobit(data_in, ldata_out,
start_ind, temp3, R, DT, MaxIters)
glm.cubic = Cubic$glm.cubic
R = Cubic$R
AIC3 = Cubic$AIC3
scaledC = Cubic$scaled
ind = which(c(AIC1, AIC2, AIC3) == min(c(AIC1, AIC2,
AIC3)))
if (ind == 1) {
mod = glm.linear
LinPart = seq(1, n2)
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
k = R[start_ind, 21]
R[start_ind, 53] = R[start_ind, 22]
xx <- sprintf("TOBIT: Fitting a linear model since it has the lowest AIC: data set %s",
code)
prof_type = 1
R[start_ind, 44] = R[start_ind, 23]
}
else if (ind == 2) {
mod = glm.quad
mu = mean(data_in)
sd = sd(data_in)
a1 = coef(mod)[[1]]
b1 = coef(mod)[[2]]
c1 = coef(mod)[[3]]
c = a1 - b1 * mu/sd + c1 * mu^2/sd^2
b = b1/sd - 2 * mu * c1/sd^2
a = c1/sd^2
if (scaledQ == 1) {
c = a1 - b1 * mu/sd + c1 * mu^2/sd^2
b = b1/sd - 2 * mu * c1/sd^2
a = c1/sd^2
}
else {
c = coef(mod)[[1]]
b = coef(mod)[[2]]
a = coef(mod)[[3]]
}
mod.fitted <- fitted.values(mod)
slopes = diff(mod.fitted)/diff(data_in)
c_points = c(data_in[1], data_in[length(data_in) -
1])
c_devs = 2 * a * c_points + b
ind_max = 2
fastest_rate = c_devs[ind_max]
ratios = fastest_rate/slopes
LinPart = c(which((ratios <= 5) & (ratios >= 0)),
n2)
ldattemp = mod.fitted[LinPart]
if (a < 0) {
xx <- sprintf("TOBIT: Fitting a quadratic (convex) model since it has the lowest AIC: data set %s",
code)
prof_type = 3
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
ldattemp = mod.fitted[LinPart]
tlag = data_in[LinPart[1]]
dattemp = data_in[LinPart]
mod <- glm(ldattemp ~ dattemp)
k = mod$coeff[[2]]
alin = mod$coeff[[1]]
}
else if (a >= 0) {
xx <- sprintf("TOBIT: Fitting a linear model since we have identified a concave quadratic (quadratic has the lowest AIC): data set %s",
code)
prof_type = 2
LinPart = seq(1, n2)
ldattemp = ldata_out[LinPart]
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
tlag = data_in[LinPart[1]]
mod <- glm.linear
k = R[start_ind, 21]
}
if (tlag == data_in[1]) {
LinPart = seq(1, n2)
LinPart = LinPart[1:(length(LinPart) - 1)]
mod <- glm.linear
k = R[start_ind, 21]
to_include_R2 = LinPart
}
else {
to_include_R2 = seq(1, length(LinPart))
}
temp3 = indexes[temp2]
R[start_ind, 44] = 1 - sum((fitted(mod)[to_include_R2] -
ldata_out[LinPart])^2)/sum((ldata_out[LinPart] -
mean(ldata_out[LinPart]))^2)
SE_temp = R[start_ind, 22]
if (tlag > data_in[1]) {
SE_temp = sqrt(diag(vcov(mod)))[[2]]
}
R[start_ind, 53] = SE_temp
R[start_ind, 50] <- mod.fitted[which(data_in == tlag)[1]] -
mod.fitted[1]
R[start_ind, 50] <- mod.fitted[which(data_in == tlag)[1]] -
mod.fitted[1]
glm.quad.fitted <- fitted.values(glm.quad)
R[start_ind, 51] <- max(glm.quad.fitted) - glm.quad.fitted[which(data_in ==
tlag)[1]]
}
else if (ind == 3) {
mod = glm.cubic
a1 = coef(mod)[[1]]
b1 = coef(mod)[[2]]
c1 = coef(mod)[[3]]
d1 = coef(mod)[[4]]
mu = mean(data_in)
sd = sd(data_in)
if (scaledC == 1) {
d = a1 - b1 * mu/sd + c1 * mu^2/sd^2 - d1 * mu^3/sd^3
c = b1/sd - 2 * mu * c1/sd^2 + 3 * mu^2 * d1/sd^3
b = c1/sd^2 - 3 * mu * d1/sd^3
a = d1/sd^3
}
else {
d = coef(mod)[[1]]
c = coef(mod)[[2]]
b = coef(mod)[[3]]
a = coef(mod)[[4]]
}
PoI = -b/(3 * a)
c_points = c(data_in[1], PoI, data_in[length(data_in) -
1])
c_devs = 3 * a * c_points^2 + 2 * b * c_points +
c
ind_min = which(abs(c_devs) == min(abs(c_devs)))
slowest_rate = c_devs[ind_min]
ind_max = which(abs(c_devs) == max(abs(c_devs)))
fastest_rate = c_devs[which(c_devs == min(c_devs))]
ind_max_neg = which(c_devs == min(c_devs))
fastest_rate = c_devs[ind_max_neg]
mod.fitted <- fitted.values(mod)
slopes = diff(mod.fitted)/diff(data_in)
ratios = fastest_rate/slopes
FastParts = which((ratios <= 5) & (ratios >= 0))
if (((PoI < data_in[1]) | (PoI > data_in[length(data_in) -
1])) & (b < 0)) {
xx <- sprintf("TOBIT: Fitting a cubic model for cubic case (d) where b <0: data set %s",
code)
fastest_rate = min(c(c_devs[1], c_devs[3]))
ratios = fastest_rate/slopes
prof_type = 7
LinPart = which((ratios <= 5) & (ratios >= 0))
temp = seq(1, n2)
LinPart = c(LinPart, temp[LinPart[length(LinPart)] +
1])
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
ldattemp = mod.fitted[LinPart]
tlag = data_in[LinPart[1]]
dattemp = data_in[LinPart]
mod <- glm(ldattemp ~ dattemp)
k = mod$coeff[[2]]
alin = mod$coeff[[1]]
}
else if (((PoI < data_in[1]) | (PoI > data_in[length(data_in) -
1])) & (b > 0)) {
xx <- sprintf("TOBIT: Fitting a linear model for cubic case (d) where b >0: data set %s",
code)
prof_type = 6
LinPart = seq(1, n2)
ldattemp = ldata_out[LinPart]
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
tlag = data_in[LinPart[1]]
mod <- glm.linear
k = R[start_ind, 21]
}
else if (ind_max_neg == 2) {
xx <- sprintf("TOBIT: Fitting a cubic model for cubic case (a): data set %s",
code)
prof_type = 4
LinPart = which((ratios <= 5) & (ratios >= 0) &
(slopes < 0))
temp = seq(1, n2)
LinPart = c(LinPart, temp[LinPart[length(LinPart)] +
1])
ldattemp = mod.fitted[LinPart]
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
ldattemp = mod.fitted[LinPart]
tlag = data_in[LinPart[1]]
dattemp = data_in[LinPart]
mod <- glm(ldattemp ~ dattemp)
k = mod$coeff[[2]]
alin = mod$coeff[[1]]
}
else if ((ind_max_neg != 2) & (FastParts[1] == 1)) {
xx <- sprintf("TOBIT: Fitting a linear model for cubic case (biii) where there is no lag phase: data set %s",
code)
prof_type = 9
LinPart = seq(1, n2)
ldattemp = ldata_out[LinPart]
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
tlag = data_in[LinPart[1]]
mod <- glm.linear
k = R[start_ind, 21]
}
else if ((ind_max_neg != 2) & (FastParts[1] != 1)) {
xx <- sprintf("TOBIT: Fitting a linear model for cubic case (bii) where there is a lag phase: data set %s",
code)
prof_type = 8
LinPart = which((ratios <= 5) & (ratios >= 0) &
(data_in[2:n2] > PoI))
temp = seq(1, n2)
LinPart = c(LinPart, temp[LinPart[length(LinPart)] +
1])
ldattemp = mod.fitted[LinPart]
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
ldattemp = mod.fitted[LinPart]
tlag = data_in[LinPart[1]]
dattemp = data_in[LinPart]
mod <- glm(ldattemp ~ dattemp)
k = mod$coeff[[2]]
alin = mod$coeff[[1]]
}
if (tlag == data_in[1]) {
LinPart = seq(1, n2)
LinPart = LinPart[1:(length(LinPart) - 1)]
mod <- glm.linear
k = R[start_ind, 21]
to_include_R2 = LinPart
}
else {
to_include_R2 = seq(1, length(LinPart))
}
temp3 = indexes[temp2]
R[start_ind, 44] = 1 - sum((fitted(mod)[to_include_R2] -
ldata_out[LinPart])^2)/sum((ldata_out[LinPart] -
mean(ldata_out[LinPart]))^2)
SE_temp = R[start_ind, 22]
if (tlag > data_in[1]) {
SE_temp = sqrt(diag(vcov(mod)))[[2]]
}
R[start_ind, 53] = SE_temp
R[start_ind, 50] <- mod.fitted[which(data_in == tlag)[1]] -
mod.fitted[1]
glm.cubic.fitted <- fitted.values(glm.cubic)
R[start_ind, 51] <- max(glm.cubic.fitted) - glm.cubic.fitted[which(data_in ==
tlag)[1]]
R[start_ind, 51] <- max(glm.cubic.fitted) - glm.cubic.fitted[which(data_in ==
tlag)[1]]
}
}
if ((plotting == TRUE) & all(ldata_out != -Inf)) {
if (imageType == "png") {
Cairo(file = filename2, type = "png")
}
else if (imageType == "pdf") {
pdf(file = filename2)
}
filename <- sprintf("patient id: %s", code)
par(mar = c(5.1, 5.1, 4.1, 2.1))
x0 = data_in[1]
xend = data_in[length(data_in)] * 1.05
y0 = 0
yend = max(ldata_out) * 1.05
if (yend == -Inf) {
yend = 10
}
xlim = c(x0, xend)
ylim = c(y0, yend)
plot(data_in, ldata_out, main = filename, xlab = "Time (hours)",
ylab = expression(log[e](parasitaemia)), cex = 2,
cex.main = 2, cex.lab = 2, cex.axis = 2, xlim = xlim,
ylim = ylim, type = "p", pch = 19)
if (R[start_ind, 56] == 0) {
if (length(fitted(mod)) == length(LinPart)) {
lines(data_in[LinPart], fitted(mod), lty = 1,
col = "blue", lwd = 3)
}
else if (prof_type == 5) {
lines(data_in[LinPart + 1], fitted(mod)[LinPart],
lty = 1, col = "blue", lwd = 3)
}
else {
lines(data_in[LinPart], fitted(mod)[LinPart],
lty = 1, col = "blue", lwd = 3)
}
}
out_temp = floor(as.numeric(R[start_ind:end_ind, 4]))
outliers = which(out_temp == 2)
if (length(outliers) > 0) {
points(data_in_store[outliers], log(data_out_store[outliers]),
col = "red", pch = 19, cex = 2)
}
if ((tlag > 0) & (R[start_ind, 56] == 0)) {
t_temp = which(data_in < tlag)
points(data_in[t_temp], log(data_out[t_temp]), col = "gray",
pch = 19, cex = 2)
}
points(data_in_store[ind_DL], log(data_out_store[ind_DL]),
col = "green", pch = 19, cex = 2)
dev.off()
}
print(xx)
if (is.null(k) == FALSE) {
mod_lin <- lm(ldata_out[-ind_DL] ~ data_in[-ind_DL])
R[start_ind, 70] = 1 - sum(residuals(mod_lin)^2)/sum((ldata_out[-ind_DL] -
mean(ldata_out[-ind_DL]))^2)
if (MAXREG == 1) {
MSRE = sqrt(mean((fitted(mod)[LinPart] - ldata_out[LinPart +
1])^2))
}
else if (tlag == data_in[1]) {
MSRE = (mean(residuals(mod)^2))^(1/2)
}
else {
MSRE = sqrt(mean((fitted(mod) - ldata_out[LinPart])^2))
}
R[start_ind, 55] = MSRE
R[temp3[LinPart], 52] = 1
if (MAXREG == 1) {
R[temp3[LinPart + 1], 45] <- fitted.values(mod)[LinPart]
R[temp3[LinPart + 1], 46] <- residuals(mod)[LinPart]
}
else if (length(fitted(mod)) == length(LinPart)) {
R[temp3[LinPart], 45] <- fitted.values(mod)
R[temp3[LinPart], 46] <- residuals(mod)
}
else {
R[temp3[LinPart], 45] <- fitted.values(mod)[LinPart]
R[temp3[LinPart], 46] <- residuals(mod)[LinPart]
}
R[start_ind, 49] <- tlag - data_in[1]
R[start_ind, 42] <- prof_type
R[start_ind, 43] <- k
knum = as.numeric(k)
int = alin + knum * as.numeric(R[start_ind, 49])
R[start_ind, 67] = int
PC50 = (log(data_out[1] * (1 - 0.5)) - alin)/knum
PC90 = (log(data_out[1] * (1 - 0.9)) - alin)/knum
PC95 = (log(data_out[1] * (1 - 0.95)) - alin)/knum
PC99 = (log(data_out[1] * (1 - 0.99)) - alin)/knum
if (PC50 > 0) {
R[start_ind, 59] = PC50
}
if (PC90 > 0) {
R[start_ind, 60] = PC90
}
if (PC95 > 0) {
R[start_ind, 61] = PC95
}
if (PC99 > 0) {
R[start_ind, 62] = PC99
}
lag_phase_attempted = 1
if ((NE == 0) | (n2 == 3)) {
lag_phase_attempted = 0
}
R[start_ind, 50] <- data_out[which(data_in == tlag)[1]] -
max(data_out[which(data_in <= tlag)])
R[start_ind, 51] <- data_out[which(data_in == tlag)[1]] -
data_out[1]
if (k > 0) {
xx <- sprintf("TOBIT: NOTE: Model is fitted, but parasites not cleared: data set %s",
code)
print(xx)
}
}
R[start_ind, 57] <- lag_phase_attempted
return(list(R = R))
}
<bytecode: 0x7816d60>
<environment: namespace:bhrcr>
--- function search by body ---
Function lagReg_tobit in namespace bhrcr has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(x_CI) == "try-error") { : the condition has length > 1
Calls: calculatePCE
Execution halted
Flavor: r-devel-linux-x86_64-debian-clang
Version: 1.0.3
Check: for non-standard things in the check directory
Result: NOTE
Found the following files/directories:
'PceEstimates'
Flavors: r-devel-linux-x86_64-debian-clang, r-devel-linux-x86_64-debian-gcc, r-devel-linux-x86_64-fedora-clang, r-devel-linux-x86_64-fedora-gcc
Version: 1.0.3
Check: examples
Result: ERROR
Running examples in ‘bhrcr-Ex.R’ failed
The error most likely occurred in:
> base::assign(".ptime", proc.time(), pos = "CheckExEnv")
> ### Name: calculatePCE
> ### Title: WWARN Parasite Clearance Estimator (PCE)
> ### Aliases: calculatePCE
>
> ### ** Examples
>
> ## Don't show:
> data("pursat")
> data = pursat[pursat["id"] <= 80 & pursat["id"] > 70,]
> output <- calculatePCE(data = data, detect.limit = 15, outlier.detect = TRUE)
Calculating WWARN PCE Estimates...
[1] "Replacing a 0 parasite level (at data point 18) with detection level: data set 71"
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
bhrcr
--- call from context ---
lagReg_tobit(R, data_in, data_out, i, start_ind, end_ind, code,
Condition_on_24_hours, fact, Threshold2, Threshold3, Detec_limit,
DT, MaxIters, plotting, Name, ind_DL, FirstHours, MaxDiffInFirstHours,
NeededInFirstHours, TimeRes1, TimeRes2, Threshold4, imageType)
--- call from argument ---
if (class(x_CI) == "try-error") {
x_CI = array(NA, c(2, 2))
} else if (is.na(x_CI[2, 1])) {
x_CI[2, 1] = NA
}
--- R stacktrace ---
where 1: lagReg_tobit(R, data_in, data_out, i, start_ind, end_ind, code,
Condition_on_24_hours, fact, Threshold2, Threshold3, Detec_limit,
DT, MaxIters, plotting, Name, ind_DL, FirstHours, MaxDiffInFirstHours,
NeededInFirstHours, TimeRes1, TimeRes2, Threshold4, imageType)
where 2: calculateParasiteClearance(Name = "pceestimates", minpara = detect.limit,
data1 = data)
where 3: withCallingHandlers(expr, message = function(c) invokeRestart("muffleMessage"))
where 4: suppressMessages(calculateParasiteClearance(Name = "pceestimates",
minpara = detect.limit, data1 = data))
where 5: withCallingHandlers(expr, warning = function(w) invokeRestart("muffleWarning"))
where 6: suppressWarnings(suppressMessages(calculateParasiteClearance(Name = "pceestimates",
minpara = detect.limit, data1 = data)))
where 7: calculatePCE(data = data, detect.limit = 15, outlier.detect = TRUE)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (R, data_in, data_out, i, start_ind, end_ind, code,
Condition_on_24_hours, fact, Threshold2, Threshold3, Detect_limit,
DT, MaxIters, plotting, Name, ind_DL, FirstHours, MaxDiffInFirstHours,
NeededInFirstHours, TimeRes1, TimeRes2, Threshold4, imageType = c("pdf",
"png"))
{
par(ask = FALSE)
imageType <- match.arg(imageType)
indexes = start_ind:end_ind
temp2 = which(R[start_ind:end_ind, 4] == 0)
temp3 = indexes[temp2]
data_out_store = data_out
data_in_store = data_in
data_in = data_in[temp2]
data_out = data_out[temp2]
ldata_out = log(data_out)
n2 = length(ldata_out)
R[temp3, 5] = ldata_out
lag_phase_attempted = 2
MAXREG = 0
filename2 <- sprintf("%s_%s.%s", Name, code, imageType)
NUM = data_in[which(data_in <= FirstHours)]
if (n2 == 3) {
NE = 0
}
else if (length(NUM) <= 1) {
NE = 1
}
else if ((max(diff(NUM)) > MaxDiffInFirstHours) | (length(NUM) <
NeededInFirstHours)) {
NE = 1
}
else {
NE = 0
}
LP_temp = which(data_out != 0)
LLP = LP_temp[length(LP_temp)]
LastPos = data_out[LLP]
k = NULL
tlag = data_in[1]
if (n2 > 3) {
Y = log(data_out[1:(n2 - 1)])
X = data_in[1:(n2 - 1)]
cf = coef(lm(Y ~ X))
x_CI = try(confint(nls(Y ~ beta * (X - x0) - log(DT),
start = c(beta = cf[[2]], x0 = -cf[[1]]/cf[[2]]))))
if (class(x_CI) == "try-error") {
x_CI = array(NA, c(2, 2))
}
else if (is.na(x_CI[2, 1])) {
x_CI[2, 1] = NA
}
}
prof_type = -9999
if (n2 < 3) {
xx <- sprintf("TOBIT: There are not enough time points - no model fitted: data set %s",
code)
R[start_ind, 56] = 1
}
else if (data_out[1] < Threshold2) {
xx <- sprintf("TOBIT: Baseline level is too low - no model fitted: data set %s",
code)
R[start_ind, 56] = 2
}
else if (n2 == 3) {
ldata_out = log(data_out)
R[temp3, 5] = ldata_out
Linear = fitLinearAndStorage(data_in, ldata_out, start_ind,
temp3, R)
glm.linear = Linear$glm.linear
R = Linear$R
AIC1 = Linear$AIC1
mod = glm.linear
LinPart = seq(1, n2)
k = R[start_ind, 21]
alin = mod$coeff[[1]]
R[start_ind, 53] = R[start_ind, 22]
xx <- sprintf("TOBIT: Fitting a linear model since it has the lowest AIC n = 3: data set %s",
code)
prof_type = 1
R[start_ind, 44] = R[start_ind, 23]
}
else if ((data_in[n2] > x_CI[[2, 1]]) & (ldata_out[n2 - 1] >
Threshold3)) {
xx <- sprintf("TOBIT: zero is not informative- no model is fitted: data set %s",
code)
R[start_ind, 56] = 3.1
}
else if (NE == 1) {
ldata_out = log(data_out)
R[temp3, 5] = ldata_out
Linear = try(fitLinearAndStorage_tobit(data_in, ldata_out,
start_ind, temp3, R, DT, MaxIters))
glm.linear = Linear$glm.linear
R = Linear$R
AIC1 = Linear$AIC1
alin = Linear$a1
mod = glm.linear
LinPart = seq(1, n2)
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
k = R[start_ind, 21]
R[start_ind, 53] = R[start_ind, 22]
xx <- sprintf("TOBIT: Fitting a linear model since it has the lowest AIC n = 4: data set %s",
code)
prof_type = 1
R[start_ind, 44] = R[start_ind, 23]
}
else if (n2 == 4) {
ldata_out = log(data_out)
R[temp3, 5] = ldata_out
Linear = try(fitLinearAndStorage_tobit(data_in, ldata_out,
start_ind, temp3, R, DT, MaxIters))
glm.linear = Linear$glm.linear
R = Linear$R
AIC1 = Linear$AIC1
alin = Linear$a1
mod = glm.linear
LinPart = seq(1, n2)
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
k = R[start_ind, 21]
R[start_ind, 53] = R[start_ind, 22]
xx <- sprintf("TOBIT: Fitting a linear model since it has the lowest AIC n = 4: data set %s",
code)
prof_type = 1
R[start_ind, 44] = R[start_ind, 23]
}
else if (n2 == 5) {
Linear = try(fitLinearAndStorage_tobit(data_in, ldata_out,
start_ind, temp3, R, DT, MaxIters))
glm.linear = Linear$glm.linear
R = Linear$R
AIC1 = Linear$AIC1
alin = Linear$a1
Quad = try(fitQuadraticAndStorage_tobit(data_in, ldata_out,
start_ind, temp3, R, DT, MaxIters))
glm.quad = Quad$glm.quad
R = Quad$R
AIC2 = Quad$AIC2
scaledQ = Quad$scaled
ind = which(c(AIC1, AIC2) == min(c(AIC1, AIC2)))
ind = which(c(AIC1, AIC2) == min(c(AIC1, AIC2)))
if (data_out[2] > (1 + fact) * data_out[1]) {
data_in_mod = data_in[2:n2]
ldata_out_mod = ldata_out[2:n2]
tlag = data_in_mod[1]
m1 <- lm(ldata_out_mod ~ poly(data_in_mod, 1, raw = TRUE))
MaxReg <- try(tobit(ldata_out_mod ~ poly(data_in_mod,
1, raw = TRUE), left = log(DT), iter.max = MaxIters))
if (class(MaxReg)[1] == "try-error") {
yy <- sprintf("Fixing error in tobit fitting, data set: %s",
code)
print(yy)
options(warn = 1)
MaxReg <- tobit(ldata_out_mod ~ poly(data_in_mod,
1, raw = TRUE), left = log(DT), init = coef(m1),
iter.max = MaxIters)
options(warn = 2)
}
AICMaxReg = AIC(MaxReg)
ss_1 = sum(residuals(glm.linear)[2:(n2 - 1)]^2)
ss_2 = sum(residuals(glm.quad)[2:(n2 - 1)]^2)
ss_3 = sum(residuals(MaxReg)[1:3]^2)
ind = which(c(ss_1, ss_2, ss_3) == min(c(ss_1, ss_2,
ss_3)))
}
if (ind == 1) {
mod = glm.linear
LinPart = seq(1, n2)
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
tlag = data_in[1]
k = R[start_ind, 21]
R[start_ind, 53] = R[start_ind, 22]
xx <- sprintf("TOBIT: Fitting a linear model since it has the lowest SS n = 5: data set %s",
code)
prof_type = 1
R[start_ind, 44] = R[start_ind, 23]
}
else if (ind == 2) {
mod = glm.quad
mu = mean(data_in)
sd = sd(data_in)
a1 = coef(mod)[[1]]
b1 = coef(mod)[[2]]
c1 = coef(mod)[[3]]
if (scaledQ == 1) {
c = a1 - b1 * mu/sd + c1 * mu^2/sd^2
b = b1/sd - 2 * mu * c1/sd^2
a = c1/sd^2
}
else {
c = coef(mod)[[1]]
b = coef(mod)[[2]]
a = coef(mod)[[3]]
}
mod.fitted <- fitted.values(mod)
slopes = diff(mod.fitted)/diff(data_in)
c_points = c(data_in[1], data_in[length(data_in) -
1])
c_devs = 2 * a * c_points + b
ind_max = 2
fastest_rate = c_devs[ind_max]
ratios = fastest_rate/slopes
LinPart = c(which((ratios <= 5) & (ratios >= 0)),
n2)
ldattemp = mod.fitted[LinPart]
if (a < 0) {
xx <- sprintf("TOBIT: Fitting a quadratic (convex) model since it has the lowest SS and n = 5: data set %s",
code)
prof_type = 3
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
ldattemp = mod.fitted[LinPart]
tlag = data_in[LinPart[1]]
dattemp = data_in[LinPart]
mod <- glm(ldattemp ~ dattemp)
k = mod$coeff[[2]]
alin = mod$coeff[[1]]
}
else if (a >= 0) {
xx <- sprintf("TOBIT: Fitting a linear model since we have identified a concave quadratic (quadratic has the lowest SS) and n = 5: data set %s",
code)
prof_type = 2
LinPart = seq(1, n2)
ldattemp = ldata_out[LinPart]
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
tlag = data_in[LinPart[1]]
mod <- glm.linear
k = R[start_ind, 21]
}
if (tlag == data_in[1]) {
LinPart = seq(1, n2)
LinPart = LinPart[1:(length(LinPart) - 1)]
mod <- glm.linear
k = R[start_ind, 21]
to_include_R2 = LinPart
}
else {
to_include_R2 = seq(1, length(LinPart))
}
temp3 = indexes[temp2]
R[start_ind, 44] = 1 - sum((fitted(mod)[to_include_R2] -
ldata_out[LinPart])^2)/sum((ldata_out[LinPart] -
mean(ldata_out[LinPart]))^2)
SE_temp = R[start_ind, 22]
if (tlag > data_in[1]) {
SE_temp = sqrt(diag(vcov(mod)))[[2]]
}
R[start_ind, 53] = SE_temp
R[start_ind, 50] <- mod.fitted[which(data_in == tlag)[1]] -
mod.fitted[1]
R[start_ind, 50] <- mod.fitted[which(data_in == tlag)[1]] -
mod.fitted[1]
glm.quad.fitted <- fitted.values(glm.quad)
R[start_ind, 51] <- max(glm.quad.fitted) - glm.quad.fitted[which(data_in ==
tlag)[1]]
}
else if (ind == 3) {
mod <- MaxReg
AIC1 = AIC(mod)
LinPart = 1:(n2 - 1)
if (LinPart[length(LinPart)] == (n2 - 1)) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
k = coef(mod)[[2]]
alin = mod$coeff[[1]]
SEs2 = sqrt(diag(vcov(mod)))
R[start_ind, 53] = SEs2[[2]]
xx <- sprintf("TOBIT: Fitting a linear model from max value since we only have 5 data points and an initial increase in parasite level (also - lowest SS): data set %s",
code)
prof_type = 5
R[start_ind, 44] = 1 - sum(residuals(mod)^2)/sum((ldata_out_mod -
mean(ldata_out_mod))^2)
MAXREG = 1
}
}
else {
Linear = fitLinearAndStorage_tobit(data_in, ldata_out,
start_ind, temp3, R, DT, MaxIters)
glm.linear = Linear$glm.linear
R = Linear$R
AIC1 = Linear$AIC1
alin = Linear$a1
Quad = fitQuadraticAndStorage_tobit(data_in, ldata_out,
start_ind, temp3, R, DT, MaxIters)
glm.quad = Quad$glm.quad
R = Quad$R
AIC2 = Quad$AIC2
scaledQ = Quad$scaled
Cubic = fitCubicAndStorage_tobit(data_in, ldata_out,
start_ind, temp3, R, DT, MaxIters)
glm.cubic = Cubic$glm.cubic
R = Cubic$R
AIC3 = Cubic$AIC3
scaledC = Cubic$scaled
ind = which(c(AIC1, AIC2, AIC3) == min(c(AIC1, AIC2,
AIC3)))
if (ind == 1) {
mod = glm.linear
LinPart = seq(1, n2)
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
k = R[start_ind, 21]
R[start_ind, 53] = R[start_ind, 22]
xx <- sprintf("TOBIT: Fitting a linear model since it has the lowest AIC: data set %s",
code)
prof_type = 1
R[start_ind, 44] = R[start_ind, 23]
}
else if (ind == 2) {
mod = glm.quad
mu = mean(data_in)
sd = sd(data_in)
a1 = coef(mod)[[1]]
b1 = coef(mod)[[2]]
c1 = coef(mod)[[3]]
c = a1 - b1 * mu/sd + c1 * mu^2/sd^2
b = b1/sd - 2 * mu * c1/sd^2
a = c1/sd^2
if (scaledQ == 1) {
c = a1 - b1 * mu/sd + c1 * mu^2/sd^2
b = b1/sd - 2 * mu * c1/sd^2
a = c1/sd^2
}
else {
c = coef(mod)[[1]]
b = coef(mod)[[2]]
a = coef(mod)[[3]]
}
mod.fitted <- fitted.values(mod)
slopes = diff(mod.fitted)/diff(data_in)
c_points = c(data_in[1], data_in[length(data_in) -
1])
c_devs = 2 * a * c_points + b
ind_max = 2
fastest_rate = c_devs[ind_max]
ratios = fastest_rate/slopes
LinPart = c(which((ratios <= 5) & (ratios >= 0)),
n2)
ldattemp = mod.fitted[LinPart]
if (a < 0) {
xx <- sprintf("TOBIT: Fitting a quadratic (convex) model since it has the lowest AIC: data set %s",
code)
prof_type = 3
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
ldattemp = mod.fitted[LinPart]
tlag = data_in[LinPart[1]]
dattemp = data_in[LinPart]
mod <- glm(ldattemp ~ dattemp)
k = mod$coeff[[2]]
alin = mod$coeff[[1]]
}
else if (a >= 0) {
xx <- sprintf("TOBIT: Fitting a linear model since we have identified a concave quadratic (quadratic has the lowest AIC): data set %s",
code)
prof_type = 2
LinPart = seq(1, n2)
ldattemp = ldata_out[LinPart]
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
tlag = data_in[LinPart[1]]
mod <- glm.linear
k = R[start_ind, 21]
}
if (tlag == data_in[1]) {
LinPart = seq(1, n2)
LinPart = LinPart[1:(length(LinPart) - 1)]
mod <- glm.linear
k = R[start_ind, 21]
to_include_R2 = LinPart
}
else {
to_include_R2 = seq(1, length(LinPart))
}
temp3 = indexes[temp2]
R[start_ind, 44] = 1 - sum((fitted(mod)[to_include_R2] -
ldata_out[LinPart])^2)/sum((ldata_out[LinPart] -
mean(ldata_out[LinPart]))^2)
SE_temp = R[start_ind, 22]
if (tlag > data_in[1]) {
SE_temp = sqrt(diag(vcov(mod)))[[2]]
}
R[start_ind, 53] = SE_temp
R[start_ind, 50] <- mod.fitted[which(data_in == tlag)[1]] -
mod.fitted[1]
R[start_ind, 50] <- mod.fitted[which(data_in == tlag)[1]] -
mod.fitted[1]
glm.quad.fitted <- fitted.values(glm.quad)
R[start_ind, 51] <- max(glm.quad.fitted) - glm.quad.fitted[which(data_in ==
tlag)[1]]
}
else if (ind == 3) {
mod = glm.cubic
a1 = coef(mod)[[1]]
b1 = coef(mod)[[2]]
c1 = coef(mod)[[3]]
d1 = coef(mod)[[4]]
mu = mean(data_in)
sd = sd(data_in)
if (scaledC == 1) {
d = a1 - b1 * mu/sd + c1 * mu^2/sd^2 - d1 * mu^3/sd^3
c = b1/sd - 2 * mu * c1/sd^2 + 3 * mu^2 * d1/sd^3
b = c1/sd^2 - 3 * mu * d1/sd^3
a = d1/sd^3
}
else {
d = coef(mod)[[1]]
c = coef(mod)[[2]]
b = coef(mod)[[3]]
a = coef(mod)[[4]]
}
PoI = -b/(3 * a)
c_points = c(data_in[1], PoI, data_in[length(data_in) -
1])
c_devs = 3 * a * c_points^2 + 2 * b * c_points +
c
ind_min = which(abs(c_devs) == min(abs(c_devs)))
slowest_rate = c_devs[ind_min]
ind_max = which(abs(c_devs) == max(abs(c_devs)))
fastest_rate = c_devs[which(c_devs == min(c_devs))]
ind_max_neg = which(c_devs == min(c_devs))
fastest_rate = c_devs[ind_max_neg]
mod.fitted <- fitted.values(mod)
slopes = diff(mod.fitted)/diff(data_in)
ratios = fastest_rate/slopes
FastParts = which((ratios <= 5) & (ratios >= 0))
if (((PoI < data_in[1]) | (PoI > data_in[length(data_in) -
1])) & (b < 0)) {
xx <- sprintf("TOBIT: Fitting a cubic model for cubic case (d) where b <0: data set %s",
code)
fastest_rate = min(c(c_devs[1], c_devs[3]))
ratios = fastest_rate/slopes
prof_type = 7
LinPart = which((ratios <= 5) & (ratios >= 0))
temp = seq(1, n2)
LinPart = c(LinPart, temp[LinPart[length(LinPart)] +
1])
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
ldattemp = mod.fitted[LinPart]
tlag = data_in[LinPart[1]]
dattemp = data_in[LinPart]
mod <- glm(ldattemp ~ dattemp)
k = mod$coeff[[2]]
alin = mod$coeff[[1]]
}
else if (((PoI < data_in[1]) | (PoI > data_in[length(data_in) -
1])) & (b > 0)) {
xx <- sprintf("TOBIT: Fitting a linear model for cubic case (d) where b >0: data set %s",
code)
prof_type = 6
LinPart = seq(1, n2)
ldattemp = ldata_out[LinPart]
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
tlag = data_in[LinPart[1]]
mod <- glm.linear
k = R[start_ind, 21]
}
else if (ind_max_neg == 2) {
xx <- sprintf("TOBIT: Fitting a cubic model for cubic case (a): data set %s",
code)
prof_type = 4
LinPart = which((ratios <= 5) & (ratios >= 0) &
(slopes < 0))
temp = seq(1, n2)
LinPart = c(LinPart, temp[LinPart[length(LinPart)] +
1])
ldattemp = mod.fitted[LinPart]
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
ldattemp = mod.fitted[LinPart]
tlag = data_in[LinPart[1]]
dattemp = data_in[LinPart]
mod <- glm(ldattemp ~ dattemp)
k = mod$coeff[[2]]
alin = mod$coeff[[1]]
}
else if ((ind_max_neg != 2) & (FastParts[1] == 1)) {
xx <- sprintf("TOBIT: Fitting a linear model for cubic case (biii) where there is no lag phase: data set %s",
code)
prof_type = 9
LinPart = seq(1, n2)
ldattemp = ldata_out[LinPart]
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
tlag = data_in[LinPart[1]]
mod <- glm.linear
k = R[start_ind, 21]
}
else if ((ind_max_neg != 2) & (FastParts[1] != 1)) {
xx <- sprintf("TOBIT: Fitting a linear model for cubic case (bii) where there is a lag phase: data set %s",
code)
prof_type = 8
LinPart = which((ratios <= 5) & (ratios >= 0) &
(data_in[2:n2] > PoI))
temp = seq(1, n2)
LinPart = c(LinPart, temp[LinPart[length(LinPart)] +
1])
ldattemp = mod.fitted[LinPart]
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
ldattemp = mod.fitted[LinPart]
tlag = data_in[LinPart[1]]
dattemp = data_in[LinPart]
mod <- glm(ldattemp ~ dattemp)
k = mod$coeff[[2]]
alin = mod$coeff[[1]]
}
if (tlag == data_in[1]) {
LinPart = seq(1, n2)
LinPart = LinPart[1:(length(LinPart) - 1)]
mod <- glm.linear
k = R[start_ind, 21]
to_include_R2 = LinPart
}
else {
to_include_R2 = seq(1, length(LinPart))
}
temp3 = indexes[temp2]
R[start_ind, 44] = 1 - sum((fitted(mod)[to_include_R2] -
ldata_out[LinPart])^2)/sum((ldata_out[LinPart] -
mean(ldata_out[LinPart]))^2)
SE_temp = R[start_ind, 22]
if (tlag > data_in[1]) {
SE_temp = sqrt(diag(vcov(mod)))[[2]]
}
R[start_ind, 53] = SE_temp
R[start_ind, 50] <- mod.fitted[which(data_in == tlag)[1]] -
mod.fitted[1]
glm.cubic.fitted <- fitted.values(glm.cubic)
R[start_ind, 51] <- max(glm.cubic.fitted) - glm.cubic.fitted[which(data_in ==
tlag)[1]]
R[start_ind, 51] <- max(glm.cubic.fitted) - glm.cubic.fitted[which(data_in ==
tlag)[1]]
}
}
if ((plotting == TRUE) & all(ldata_out != -Inf)) {
if (imageType == "png") {
Cairo(file = filename2, type = "png")
}
else if (imageType == "pdf") {
pdf(file = filename2)
}
filename <- sprintf("patient id: %s", code)
par(mar = c(5.1, 5.1, 4.1, 2.1))
x0 = data_in[1]
xend = data_in[length(data_in)] * 1.05
y0 = 0
yend = max(ldata_out) * 1.05
if (yend == -Inf) {
yend = 10
}
xlim = c(x0, xend)
ylim = c(y0, yend)
plot(data_in, ldata_out, main = filename, xlab = "Time (hours)",
ylab = expression(log[e](parasitaemia)), cex = 2,
cex.main = 2, cex.lab = 2, cex.axis = 2, xlim = xlim,
ylim = ylim, type = "p", pch = 19)
if (R[start_ind, 56] == 0) {
if (length(fitted(mod)) == length(LinPart)) {
lines(data_in[LinPart], fitted(mod), lty = 1,
col = "blue", lwd = 3)
}
else if (prof_type == 5) {
lines(data_in[LinPart + 1], fitted(mod)[LinPart],
lty = 1, col = "blue", lwd = 3)
}
else {
lines(data_in[LinPart], fitted(mod)[LinPart],
lty = 1, col = "blue", lwd = 3)
}
}
out_temp = floor(as.numeric(R[start_ind:end_ind, 4]))
outliers = which(out_temp == 2)
if (length(outliers) > 0) {
points(data_in_store[outliers], log(data_out_store[outliers]),
col = "red", pch = 19, cex = 2)
}
if ((tlag > 0) & (R[start_ind, 56] == 0)) {
t_temp = which(data_in < tlag)
points(data_in[t_temp], log(data_out[t_temp]), col = "gray",
pch = 19, cex = 2)
}
points(data_in_store[ind_DL], log(data_out_store[ind_DL]),
col = "green", pch = 19, cex = 2)
dev.off()
}
print(xx)
if (is.null(k) == FALSE) {
mod_lin <- lm(ldata_out[-ind_DL] ~ data_in[-ind_DL])
R[start_ind, 70] = 1 - sum(residuals(mod_lin)^2)/sum((ldata_out[-ind_DL] -
mean(ldata_out[-ind_DL]))^2)
if (MAXREG == 1) {
MSRE = sqrt(mean((fitted(mod)[LinPart] - ldata_out[LinPart +
1])^2))
}
else if (tlag == data_in[1]) {
MSRE = (mean(residuals(mod)^2))^(1/2)
}
else {
MSRE = sqrt(mean((fitted(mod) - ldata_out[LinPart])^2))
}
R[start_ind, 55] = MSRE
R[temp3[LinPart], 52] = 1
if (MAXREG == 1) {
R[temp3[LinPart + 1], 45] <- fitted.values(mod)[LinPart]
R[temp3[LinPart + 1], 46] <- residuals(mod)[LinPart]
}
else if (length(fitted(mod)) == length(LinPart)) {
R[temp3[LinPart], 45] <- fitted.values(mod)
R[temp3[LinPart], 46] <- residuals(mod)
}
else {
R[temp3[LinPart], 45] <- fitted.values(mod)[LinPart]
R[temp3[LinPart], 46] <- residuals(mod)[LinPart]
}
R[start_ind, 49] <- tlag - data_in[1]
R[start_ind, 42] <- prof_type
R[start_ind, 43] <- k
knum = as.numeric(k)
int = alin + knum * as.numeric(R[start_ind, 49])
R[start_ind, 67] = int
PC50 = (log(data_out[1] * (1 - 0.5)) - alin)/knum
PC90 = (log(data_out[1] * (1 - 0.9)) - alin)/knum
PC95 = (log(data_out[1] * (1 - 0.95)) - alin)/knum
PC99 = (log(data_out[1] * (1 - 0.99)) - alin)/knum
if (PC50 > 0) {
R[start_ind, 59] = PC50
}
if (PC90 > 0) {
R[start_ind, 60] = PC90
}
if (PC95 > 0) {
R[start_ind, 61] = PC95
}
if (PC99 > 0) {
R[start_ind, 62] = PC99
}
lag_phase_attempted = 1
if ((NE == 0) | (n2 == 3)) {
lag_phase_attempted = 0
}
R[start_ind, 50] <- data_out[which(data_in == tlag)[1]] -
max(data_out[which(data_in <= tlag)])
R[start_ind, 51] <- data_out[which(data_in == tlag)[1]] -
data_out[1]
if (k > 0) {
xx <- sprintf("TOBIT: NOTE: Model is fitted, but parasites not cleared: data set %s",
code)
print(xx)
}
}
R[start_ind, 57] <- lag_phase_attempted
return(list(R = R))
}
<bytecode: 0x5627aa43bb38>
<environment: namespace:bhrcr>
--- function search by body ---
Function lagReg_tobit in namespace bhrcr has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(x_CI) == "try-error") { : the condition has length > 1
Calls: calculatePCE
Execution halted
Flavor: r-devel-linux-x86_64-debian-gcc
Version: 1.0.3
Check: examples
Result: ERROR
Running examples in ‘bhrcr-Ex.R’ failed
The error most likely occurred in:
> ### Name: calculatePCE
> ### Title: WWARN Parasite Clearance Estimator (PCE)
> ### Aliases: calculatePCE
>
> ### ** Examples
>
> ## Don't show:
> data("pursat")
> data = pursat[pursat["id"] <= 80 & pursat["id"] > 70,]
> output <- calculatePCE(data = data, detect.limit = 15, outlier.detect = TRUE)
Calculating WWARN PCE Estimates...
[1] "Replacing a 0 parasite level (at data point 18) with detection level: data set 71"
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
bhrcr
--- call from context ---
lagReg_tobit(R, data_in, data_out, i, start_ind, end_ind, code,
Condition_on_24_hours, fact, Threshold2, Threshold3, Detec_limit,
DT, MaxIters, plotting, Name, ind_DL, FirstHours, MaxDiffInFirstHours,
NeededInFirstHours, TimeRes1, TimeRes2, Threshold4, imageType)
--- call from argument ---
if (class(x_CI) == "try-error") {
x_CI = array(NA, c(2, 2))
} else if (is.na(x_CI[2, 1])) {
x_CI[2, 1] = NA
}
--- R stacktrace ---
where 1: lagReg_tobit(R, data_in, data_out, i, start_ind, end_ind, code,
Condition_on_24_hours, fact, Threshold2, Threshold3, Detec_limit,
DT, MaxIters, plotting, Name, ind_DL, FirstHours, MaxDiffInFirstHours,
NeededInFirstHours, TimeRes1, TimeRes2, Threshold4, imageType)
where 2: calculateParasiteClearance(Name = "pceestimates", minpara = detect.limit,
data1 = data)
where 3: withCallingHandlers(expr, message = function(c) invokeRestart("muffleMessage"))
where 4: suppressMessages(calculateParasiteClearance(Name = "pceestimates",
minpara = detect.limit, data1 = data))
where 5: withCallingHandlers(expr, warning = function(w) invokeRestart("muffleWarning"))
where 6: suppressWarnings(suppressMessages(calculateParasiteClearance(Name = "pceestimates",
minpara = detect.limit, data1 = data)))
where 7: calculatePCE(data = data, detect.limit = 15, outlier.detect = TRUE)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (R, data_in, data_out, i, start_ind, end_ind, code,
Condition_on_24_hours, fact, Threshold2, Threshold3, Detect_limit,
DT, MaxIters, plotting, Name, ind_DL, FirstHours, MaxDiffInFirstHours,
NeededInFirstHours, TimeRes1, TimeRes2, Threshold4, imageType = c("pdf",
"png"))
{
par(ask = FALSE)
imageType <- match.arg(imageType)
indexes = start_ind:end_ind
temp2 = which(R[start_ind:end_ind, 4] == 0)
temp3 = indexes[temp2]
data_out_store = data_out
data_in_store = data_in
data_in = data_in[temp2]
data_out = data_out[temp2]
ldata_out = log(data_out)
n2 = length(ldata_out)
R[temp3, 5] = ldata_out
lag_phase_attempted = 2
MAXREG = 0
filename2 <- sprintf("%s_%s.%s", Name, code, imageType)
NUM = data_in[which(data_in <= FirstHours)]
if (n2 == 3) {
NE = 0
}
else if (length(NUM) <= 1) {
NE = 1
}
else if ((max(diff(NUM)) > MaxDiffInFirstHours) | (length(NUM) <
NeededInFirstHours)) {
NE = 1
}
else {
NE = 0
}
LP_temp = which(data_out != 0)
LLP = LP_temp[length(LP_temp)]
LastPos = data_out[LLP]
k = NULL
tlag = data_in[1]
if (n2 > 3) {
Y = log(data_out[1:(n2 - 1)])
X = data_in[1:(n2 - 1)]
cf = coef(lm(Y ~ X))
x_CI = try(confint(nls(Y ~ beta * (X - x0) - log(DT),
start = c(beta = cf[[2]], x0 = -cf[[1]]/cf[[2]]))))
if (class(x_CI) == "try-error") {
x_CI = array(NA, c(2, 2))
}
else if (is.na(x_CI[2, 1])) {
x_CI[2, 1] = NA
}
}
prof_type = -9999
if (n2 < 3) {
xx <- sprintf("TOBIT: There are not enough time points - no model fitted: data set %s",
code)
R[start_ind, 56] = 1
}
else if (data_out[1] < Threshold2) {
xx <- sprintf("TOBIT: Baseline level is too low - no model fitted: data set %s",
code)
R[start_ind, 56] = 2
}
else if (n2 == 3) {
ldata_out = log(data_out)
R[temp3, 5] = ldata_out
Linear = fitLinearAndStorage(data_in, ldata_out, start_ind,
temp3, R)
glm.linear = Linear$glm.linear
R = Linear$R
AIC1 = Linear$AIC1
mod = glm.linear
LinPart = seq(1, n2)
k = R[start_ind, 21]
alin = mod$coeff[[1]]
R[start_ind, 53] = R[start_ind, 22]
xx <- sprintf("TOBIT: Fitting a linear model since it has the lowest AIC n = 3: data set %s",
code)
prof_type = 1
R[start_ind, 44] = R[start_ind, 23]
}
else if ((data_in[n2] > x_CI[[2, 1]]) & (ldata_out[n2 - 1] >
Threshold3)) {
xx <- sprintf("TOBIT: zero is not informative- no model is fitted: data set %s",
code)
R[start_ind, 56] = 3.1
}
else if (NE == 1) {
ldata_out = log(data_out)
R[temp3, 5] = ldata_out
Linear = try(fitLinearAndStorage_tobit(data_in, ldata_out,
start_ind, temp3, R, DT, MaxIters))
glm.linear = Linear$glm.linear
R = Linear$R
AIC1 = Linear$AIC1
alin = Linear$a1
mod = glm.linear
LinPart = seq(1, n2)
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
k = R[start_ind, 21]
R[start_ind, 53] = R[start_ind, 22]
xx <- sprintf("TOBIT: Fitting a linear model since it has the lowest AIC n = 4: data set %s",
code)
prof_type = 1
R[start_ind, 44] = R[start_ind, 23]
}
else if (n2 == 4) {
ldata_out = log(data_out)
R[temp3, 5] = ldata_out
Linear = try(fitLinearAndStorage_tobit(data_in, ldata_out,
start_ind, temp3, R, DT, MaxIters))
glm.linear = Linear$glm.linear
R = Linear$R
AIC1 = Linear$AIC1
alin = Linear$a1
mod = glm.linear
LinPart = seq(1, n2)
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
k = R[start_ind, 21]
R[start_ind, 53] = R[start_ind, 22]
xx <- sprintf("TOBIT: Fitting a linear model since it has the lowest AIC n = 4: data set %s",
code)
prof_type = 1
R[start_ind, 44] = R[start_ind, 23]
}
else if (n2 == 5) {
Linear = try(fitLinearAndStorage_tobit(data_in, ldata_out,
start_ind, temp3, R, DT, MaxIters))
glm.linear = Linear$glm.linear
R = Linear$R
AIC1 = Linear$AIC1
alin = Linear$a1
Quad = try(fitQuadraticAndStorage_tobit(data_in, ldata_out,
start_ind, temp3, R, DT, MaxIters))
glm.quad = Quad$glm.quad
R = Quad$R
AIC2 = Quad$AIC2
scaledQ = Quad$scaled
ind = which(c(AIC1, AIC2) == min(c(AIC1, AIC2)))
ind = which(c(AIC1, AIC2) == min(c(AIC1, AIC2)))
if (data_out[2] > (1 + fact) * data_out[1]) {
data_in_mod = data_in[2:n2]
ldata_out_mod = ldata_out[2:n2]
tlag = data_in_mod[1]
m1 <- lm(ldata_out_mod ~ poly(data_in_mod, 1, raw = TRUE))
MaxReg <- try(tobit(ldata_out_mod ~ poly(data_in_mod,
1, raw = TRUE), left = log(DT), iter.max = MaxIters))
if (class(MaxReg)[1] == "try-error") {
yy <- sprintf("Fixing error in tobit fitting, data set: %s",
code)
print(yy)
options(warn = 1)
MaxReg <- tobit(ldata_out_mod ~ poly(data_in_mod,
1, raw = TRUE), left = log(DT), init = coef(m1),
iter.max = MaxIters)
options(warn = 2)
}
AICMaxReg = AIC(MaxReg)
ss_1 = sum(residuals(glm.linear)[2:(n2 - 1)]^2)
ss_2 = sum(residuals(glm.quad)[2:(n2 - 1)]^2)
ss_3 = sum(residuals(MaxReg)[1:3]^2)
ind = which(c(ss_1, ss_2, ss_3) == min(c(ss_1, ss_2,
ss_3)))
}
if (ind == 1) {
mod = glm.linear
LinPart = seq(1, n2)
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
tlag = data_in[1]
k = R[start_ind, 21]
R[start_ind, 53] = R[start_ind, 22]
xx <- sprintf("TOBIT: Fitting a linear model since it has the lowest SS n = 5: data set %s",
code)
prof_type = 1
R[start_ind, 44] = R[start_ind, 23]
}
else if (ind == 2) {
mod = glm.quad
mu = mean(data_in)
sd = sd(data_in)
a1 = coef(mod)[[1]]
b1 = coef(mod)[[2]]
c1 = coef(mod)[[3]]
if (scaledQ == 1) {
c = a1 - b1 * mu/sd + c1 * mu^2/sd^2
b = b1/sd - 2 * mu * c1/sd^2
a = c1/sd^2
}
else {
c = coef(mod)[[1]]
b = coef(mod)[[2]]
a = coef(mod)[[3]]
}
mod.fitted <- fitted.values(mod)
slopes = diff(mod.fitted)/diff(data_in)
c_points = c(data_in[1], data_in[length(data_in) -
1])
c_devs = 2 * a * c_points + b
ind_max = 2
fastest_rate = c_devs[ind_max]
ratios = fastest_rate/slopes
LinPart = c(which((ratios <= 5) & (ratios >= 0)),
n2)
ldattemp = mod.fitted[LinPart]
if (a < 0) {
xx <- sprintf("TOBIT: Fitting a quadratic (convex) model since it has the lowest SS and n = 5: data set %s",
code)
prof_type = 3
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
ldattemp = mod.fitted[LinPart]
tlag = data_in[LinPart[1]]
dattemp = data_in[LinPart]
mod <- glm(ldattemp ~ dattemp)
k = mod$coeff[[2]]
alin = mod$coeff[[1]]
}
else if (a >= 0) {
xx <- sprintf("TOBIT: Fitting a linear model since we have identified a concave quadratic (quadratic has the lowest SS) and n = 5: data set %s",
code)
prof_type = 2
LinPart = seq(1, n2)
ldattemp = ldata_out[LinPart]
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
tlag = data_in[LinPart[1]]
mod <- glm.linear
k = R[start_ind, 21]
}
if (tlag == data_in[1]) {
LinPart = seq(1, n2)
LinPart = LinPart[1:(length(LinPart) - 1)]
mod <- glm.linear
k = R[start_ind, 21]
to_include_R2 = LinPart
}
else {
to_include_R2 = seq(1, length(LinPart))
}
temp3 = indexes[temp2]
R[start_ind, 44] = 1 - sum((fitted(mod)[to_include_R2] -
ldata_out[LinPart])^2)/sum((ldata_out[LinPart] -
mean(ldata_out[LinPart]))^2)
SE_temp = R[start_ind, 22]
if (tlag > data_in[1]) {
SE_temp = sqrt(diag(vcov(mod)))[[2]]
}
R[start_ind, 53] = SE_temp
R[start_ind, 50] <- mod.fitted[which(data_in == tlag)[1]] -
mod.fitted[1]
R[start_ind, 50] <- mod.fitted[which(data_in == tlag)[1]] -
mod.fitted[1]
glm.quad.fitted <- fitted.values(glm.quad)
R[start_ind, 51] <- max(glm.quad.fitted) - glm.quad.fitted[which(data_in ==
tlag)[1]]
}
else if (ind == 3) {
mod <- MaxReg
AIC1 = AIC(mod)
LinPart = 1:(n2 - 1)
if (LinPart[length(LinPart)] == (n2 - 1)) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
k = coef(mod)[[2]]
alin = mod$coeff[[1]]
SEs2 = sqrt(diag(vcov(mod)))
R[start_ind, 53] = SEs2[[2]]
xx <- sprintf("TOBIT: Fitting a linear model from max value since we only have 5 data points and an initial increase in parasite level (also - lowest SS): data set %s",
code)
prof_type = 5
R[start_ind, 44] = 1 - sum(residuals(mod)^2)/sum((ldata_out_mod -
mean(ldata_out_mod))^2)
MAXREG = 1
}
}
else {
Linear = fitLinearAndStorage_tobit(data_in, ldata_out,
start_ind, temp3, R, DT, MaxIters)
glm.linear = Linear$glm.linear
R = Linear$R
AIC1 = Linear$AIC1
alin = Linear$a1
Quad = fitQuadraticAndStorage_tobit(data_in, ldata_out,
start_ind, temp3, R, DT, MaxIters)
glm.quad = Quad$glm.quad
R = Quad$R
AIC2 = Quad$AIC2
scaledQ = Quad$scaled
Cubic = fitCubicAndStorage_tobit(data_in, ldata_out,
start_ind, temp3, R, DT, MaxIters)
glm.cubic = Cubic$glm.cubic
R = Cubic$R
AIC3 = Cubic$AIC3
scaledC = Cubic$scaled
ind = which(c(AIC1, AIC2, AIC3) == min(c(AIC1, AIC2,
AIC3)))
if (ind == 1) {
mod = glm.linear
LinPart = seq(1, n2)
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
k = R[start_ind, 21]
R[start_ind, 53] = R[start_ind, 22]
xx <- sprintf("TOBIT: Fitting a linear model since it has the lowest AIC: data set %s",
code)
prof_type = 1
R[start_ind, 44] = R[start_ind, 23]
}
else if (ind == 2) {
mod = glm.quad
mu = mean(data_in)
sd = sd(data_in)
a1 = coef(mod)[[1]]
b1 = coef(mod)[[2]]
c1 = coef(mod)[[3]]
c = a1 - b1 * mu/sd + c1 * mu^2/sd^2
b = b1/sd - 2 * mu * c1/sd^2
a = c1/sd^2
if (scaledQ == 1) {
c = a1 - b1 * mu/sd + c1 * mu^2/sd^2
b = b1/sd - 2 * mu * c1/sd^2
a = c1/sd^2
}
else {
c = coef(mod)[[1]]
b = coef(mod)[[2]]
a = coef(mod)[[3]]
}
mod.fitted <- fitted.values(mod)
slopes = diff(mod.fitted)/diff(data_in)
c_points = c(data_in[1], data_in[length(data_in) -
1])
c_devs = 2 * a * c_points + b
ind_max = 2
fastest_rate = c_devs[ind_max]
ratios = fastest_rate/slopes
LinPart = c(which((ratios <= 5) & (ratios >= 0)),
n2)
ldattemp = mod.fitted[LinPart]
if (a < 0) {
xx <- sprintf("TOBIT: Fitting a quadratic (convex) model since it has the lowest AIC: data set %s",
code)
prof_type = 3
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
ldattemp = mod.fitted[LinPart]
tlag = data_in[LinPart[1]]
dattemp = data_in[LinPart]
mod <- glm(ldattemp ~ dattemp)
k = mod$coeff[[2]]
alin = mod$coeff[[1]]
}
else if (a >= 0) {
xx <- sprintf("TOBIT: Fitting a linear model since we have identified a concave quadratic (quadratic has the lowest AIC): data set %s",
code)
prof_type = 2
LinPart = seq(1, n2)
ldattemp = ldata_out[LinPart]
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
tlag = data_in[LinPart[1]]
mod <- glm.linear
k = R[start_ind, 21]
}
if (tlag == data_in[1]) {
LinPart = seq(1, n2)
LinPart = LinPart[1:(length(LinPart) - 1)]
mod <- glm.linear
k = R[start_ind, 21]
to_include_R2 = LinPart
}
else {
to_include_R2 = seq(1, length(LinPart))
}
temp3 = indexes[temp2]
R[start_ind, 44] = 1 - sum((fitted(mod)[to_include_R2] -
ldata_out[LinPart])^2)/sum((ldata_out[LinPart] -
mean(ldata_out[LinPart]))^2)
SE_temp = R[start_ind, 22]
if (tlag > data_in[1]) {
SE_temp = sqrt(diag(vcov(mod)))[[2]]
}
R[start_ind, 53] = SE_temp
R[start_ind, 50] <- mod.fitted[which(data_in == tlag)[1]] -
mod.fitted[1]
R[start_ind, 50] <- mod.fitted[which(data_in == tlag)[1]] -
mod.fitted[1]
glm.quad.fitted <- fitted.values(glm.quad)
R[start_ind, 51] <- max(glm.quad.fitted) - glm.quad.fitted[which(data_in ==
tlag)[1]]
}
else if (ind == 3) {
mod = glm.cubic
a1 = coef(mod)[[1]]
b1 = coef(mod)[[2]]
c1 = coef(mod)[[3]]
d1 = coef(mod)[[4]]
mu = mean(data_in)
sd = sd(data_in)
if (scaledC == 1) {
d = a1 - b1 * mu/sd + c1 * mu^2/sd^2 - d1 * mu^3/sd^3
c = b1/sd - 2 * mu * c1/sd^2 + 3 * mu^2 * d1/sd^3
b = c1/sd^2 - 3 * mu * d1/sd^3
a = d1/sd^3
}
else {
d = coef(mod)[[1]]
c = coef(mod)[[2]]
b = coef(mod)[[3]]
a = coef(mod)[[4]]
}
PoI = -b/(3 * a)
c_points = c(data_in[1], PoI, data_in[length(data_in) -
1])
c_devs = 3 * a * c_points^2 + 2 * b * c_points +
c
ind_min = which(abs(c_devs) == min(abs(c_devs)))
slowest_rate = c_devs[ind_min]
ind_max = which(abs(c_devs) == max(abs(c_devs)))
fastest_rate = c_devs[which(c_devs == min(c_devs))]
ind_max_neg = which(c_devs == min(c_devs))
fastest_rate = c_devs[ind_max_neg]
mod.fitted <- fitted.values(mod)
slopes = diff(mod.fitted)/diff(data_in)
ratios = fastest_rate/slopes
FastParts = which((ratios <= 5) & (ratios >= 0))
if (((PoI < data_in[1]) | (PoI > data_in[length(data_in) -
1])) & (b < 0)) {
xx <- sprintf("TOBIT: Fitting a cubic model for cubic case (d) where b <0: data set %s",
code)
fastest_rate = min(c(c_devs[1], c_devs[3]))
ratios = fastest_rate/slopes
prof_type = 7
LinPart = which((ratios <= 5) & (ratios >= 0))
temp = seq(1, n2)
LinPart = c(LinPart, temp[LinPart[length(LinPart)] +
1])
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
ldattemp = mod.fitted[LinPart]
tlag = data_in[LinPart[1]]
dattemp = data_in[LinPart]
mod <- glm(ldattemp ~ dattemp)
k = mod$coeff[[2]]
alin = mod$coeff[[1]]
}
else if (((PoI < data_in[1]) | (PoI > data_in[length(data_in) -
1])) & (b > 0)) {
xx <- sprintf("TOBIT: Fitting a linear model for cubic case (d) where b >0: data set %s",
code)
prof_type = 6
LinPart = seq(1, n2)
ldattemp = ldata_out[LinPart]
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
tlag = data_in[LinPart[1]]
mod <- glm.linear
k = R[start_ind, 21]
}
else if (ind_max_neg == 2) {
xx <- sprintf("TOBIT: Fitting a cubic model for cubic case (a): data set %s",
code)
prof_type = 4
LinPart = which((ratios <= 5) & (ratios >= 0) &
(slopes < 0))
temp = seq(1, n2)
LinPart = c(LinPart, temp[LinPart[length(LinPart)] +
1])
ldattemp = mod.fitted[LinPart]
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
ldattemp = mod.fitted[LinPart]
tlag = data_in[LinPart[1]]
dattemp = data_in[LinPart]
mod <- glm(ldattemp ~ dattemp)
k = mod$coeff[[2]]
alin = mod$coeff[[1]]
}
else if ((ind_max_neg != 2) & (FastParts[1] == 1)) {
xx <- sprintf("TOBIT: Fitting a linear model for cubic case (biii) where there is no lag phase: data set %s",
code)
prof_type = 9
LinPart = seq(1, n2)
ldattemp = ldata_out[LinPart]
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
tlag = data_in[LinPart[1]]
mod <- glm.linear
k = R[start_ind, 21]
}
else if ((ind_max_neg != 2) & (FastParts[1] != 1)) {
xx <- sprintf("TOBIT: Fitting a linear model for cubic case (bii) where there is a lag phase: data set %s",
code)
prof_type = 8
LinPart = which((ratios <= 5) & (ratios >= 0) &
(data_in[2:n2] > PoI))
temp = seq(1, n2)
LinPart = c(LinPart, temp[LinPart[length(LinPart)] +
1])
ldattemp = mod.fitted[LinPart]
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
ldattemp = mod.fitted[LinPart]
tlag = data_in[LinPart[1]]
dattemp = data_in[LinPart]
mod <- glm(ldattemp ~ dattemp)
k = mod$coeff[[2]]
alin = mod$coeff[[1]]
}
if (tlag == data_in[1]) {
LinPart = seq(1, n2)
LinPart = LinPart[1:(length(LinPart) - 1)]
mod <- glm.linear
k = R[start_ind, 21]
to_include_R2 = LinPart
}
else {
to_include_R2 = seq(1, length(LinPart))
}
temp3 = indexes[temp2]
R[start_ind, 44] = 1 - sum((fitted(mod)[to_include_R2] -
ldata_out[LinPart])^2)/sum((ldata_out[LinPart] -
mean(ldata_out[LinPart]))^2)
SE_temp = R[start_ind, 22]
if (tlag > data_in[1]) {
SE_temp = sqrt(diag(vcov(mod)))[[2]]
}
R[start_ind, 53] = SE_temp
R[start_ind, 50] <- mod.fitted[which(data_in == tlag)[1]] -
mod.fitted[1]
glm.cubic.fitted <- fitted.values(glm.cubic)
R[start_ind, 51] <- max(glm.cubic.fitted) - glm.cubic.fitted[which(data_in ==
tlag)[1]]
R[start_ind, 51] <- max(glm.cubic.fitted) - glm.cubic.fitted[which(data_in ==
tlag)[1]]
}
}
if ((plotting == TRUE) & all(ldata_out != -Inf)) {
if (imageType == "png") {
Cairo(file = filename2, type = "png")
}
else if (imageType == "pdf") {
pdf(file = filename2)
}
filename <- sprintf("patient id: %s", code)
par(mar = c(5.1, 5.1, 4.1, 2.1))
x0 = data_in[1]
xend = data_in[length(data_in)] * 1.05
y0 = 0
yend = max(ldata_out) * 1.05
if (yend == -Inf) {
yend = 10
}
xlim = c(x0, xend)
ylim = c(y0, yend)
plot(data_in, ldata_out, main = filename, xlab = "Time (hours)",
ylab = expression(log[e](parasitaemia)), cex = 2,
cex.main = 2, cex.lab = 2, cex.axis = 2, xlim = xlim,
ylim = ylim, type = "p", pch = 19)
if (R[start_ind, 56] == 0) {
if (length(fitted(mod)) == length(LinPart)) {
lines(data_in[LinPart], fitted(mod), lty = 1,
col = "blue", lwd = 3)
}
else if (prof_type == 5) {
lines(data_in[LinPart + 1], fitted(mod)[LinPart],
lty = 1, col = "blue", lwd = 3)
}
else {
lines(data_in[LinPart], fitted(mod)[LinPart],
lty = 1, col = "blue", lwd = 3)
}
}
out_temp = floor(as.numeric(R[start_ind:end_ind, 4]))
outliers = which(out_temp == 2)
if (length(outliers) > 0) {
points(data_in_store[outliers], log(data_out_store[outliers]),
col = "red", pch = 19, cex = 2)
}
if ((tlag > 0) & (R[start_ind, 56] == 0)) {
t_temp = which(data_in < tlag)
points(data_in[t_temp], log(data_out[t_temp]), col = "gray",
pch = 19, cex = 2)
}
points(data_in_store[ind_DL], log(data_out_store[ind_DL]),
col = "green", pch = 19, cex = 2)
dev.off()
}
print(xx)
if (is.null(k) == FALSE) {
mod_lin <- lm(ldata_out[-ind_DL] ~ data_in[-ind_DL])
R[start_ind, 70] = 1 - sum(residuals(mod_lin)^2)/sum((ldata_out[-ind_DL] -
mean(ldata_out[-ind_DL]))^2)
if (MAXREG == 1) {
MSRE = sqrt(mean((fitted(mod)[LinPart] - ldata_out[LinPart +
1])^2))
}
else if (tlag == data_in[1]) {
MSRE = (mean(residuals(mod)^2))^(1/2)
}
else {
MSRE = sqrt(mean((fitted(mod) - ldata_out[LinPart])^2))
}
R[start_ind, 55] = MSRE
R[temp3[LinPart], 52] = 1
if (MAXREG == 1) {
R[temp3[LinPart + 1], 45] <- fitted.values(mod)[LinPart]
R[temp3[LinPart + 1], 46] <- residuals(mod)[LinPart]
}
else if (length(fitted(mod)) == length(LinPart)) {
R[temp3[LinPart], 45] <- fitted.values(mod)
R[temp3[LinPart], 46] <- residuals(mod)
}
else {
R[temp3[LinPart], 45] <- fitted.values(mod)[LinPart]
R[temp3[LinPart], 46] <- residuals(mod)[LinPart]
}
R[start_ind, 49] <- tlag - data_in[1]
R[start_ind, 42] <- prof_type
R[start_ind, 43] <- k
knum = as.numeric(k)
int = alin + knum * as.numeric(R[start_ind, 49])
R[start_ind, 67] = int
PC50 = (log(data_out[1] * (1 - 0.5)) - alin)/knum
PC90 = (log(data_out[1] * (1 - 0.9)) - alin)/knum
PC95 = (log(data_out[1] * (1 - 0.95)) - alin)/knum
PC99 = (log(data_out[1] * (1 - 0.99)) - alin)/knum
if (PC50 > 0) {
R[start_ind, 59] = PC50
}
if (PC90 > 0) {
R[start_ind, 60] = PC90
}
if (PC95 > 0) {
R[start_ind, 61] = PC95
}
if (PC99 > 0) {
R[start_ind, 62] = PC99
}
lag_phase_attempted = 1
if ((NE == 0) | (n2 == 3)) {
lag_phase_attempted = 0
}
R[start_ind, 50] <- data_out[which(data_in == tlag)[1]] -
max(data_out[which(data_in <= tlag)])
R[start_ind, 51] <- data_out[which(data_in == tlag)[1]] -
data_out[1]
if (k > 0) {
xx <- sprintf("TOBIT: NOTE: Model is fitted, but parasites not cleared: data set %s",
code)
print(xx)
}
}
R[start_ind, 57] <- lag_phase_attempted
return(list(R = R))
}
<bytecode: 0x55b3528>
<environment: namespace:bhrcr>
--- function search by body ---
Function lagReg_tobit in namespace bhrcr has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(x_CI) == "try-error") { : the condition has length > 1
Calls: calculatePCE
Execution halted
Flavor: r-devel-linux-x86_64-fedora-clang
Version: 1.0.3
Check: examples
Result: ERROR
Running examples in ‘bhrcr-Ex.R’ failed
The error most likely occurred in:
> ### Name: calculatePCE
> ### Title: WWARN Parasite Clearance Estimator (PCE)
> ### Aliases: calculatePCE
>
> ### ** Examples
>
> ## Don't show:
> data("pursat")
> data = pursat[pursat["id"] <= 80 & pursat["id"] > 70,]
> output <- calculatePCE(data = data, detect.limit = 15, outlier.detect = TRUE)
Calculating WWARN PCE Estimates...
[1] "Replacing a 0 parasite level (at data point 18) with detection level: data set 71"
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
bhrcr
--- call from context ---
lagReg_tobit(R, data_in, data_out, i, start_ind, end_ind, code,
Condition_on_24_hours, fact, Threshold2, Threshold3, Detec_limit,
DT, MaxIters, plotting, Name, ind_DL, FirstHours, MaxDiffInFirstHours,
NeededInFirstHours, TimeRes1, TimeRes2, Threshold4, imageType)
--- call from argument ---
if (class(x_CI) == "try-error") {
x_CI = array(NA, c(2, 2))
} else if (is.na(x_CI[2, 1])) {
x_CI[2, 1] = NA
}
--- R stacktrace ---
where 1: lagReg_tobit(R, data_in, data_out, i, start_ind, end_ind, code,
Condition_on_24_hours, fact, Threshold2, Threshold3, Detec_limit,
DT, MaxIters, plotting, Name, ind_DL, FirstHours, MaxDiffInFirstHours,
NeededInFirstHours, TimeRes1, TimeRes2, Threshold4, imageType)
where 2: calculateParasiteClearance(Name = "pceestimates", minpara = detect.limit,
data1 = data)
where 3: withCallingHandlers(expr, message = function(c) invokeRestart("muffleMessage"))
where 4: suppressMessages(calculateParasiteClearance(Name = "pceestimates",
minpara = detect.limit, data1 = data))
where 5: withCallingHandlers(expr, warning = function(w) invokeRestart("muffleWarning"))
where 6: suppressWarnings(suppressMessages(calculateParasiteClearance(Name = "pceestimates",
minpara = detect.limit, data1 = data)))
where 7: calculatePCE(data = data, detect.limit = 15, outlier.detect = TRUE)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (R, data_in, data_out, i, start_ind, end_ind, code,
Condition_on_24_hours, fact, Threshold2, Threshold3, Detect_limit,
DT, MaxIters, plotting, Name, ind_DL, FirstHours, MaxDiffInFirstHours,
NeededInFirstHours, TimeRes1, TimeRes2, Threshold4, imageType = c("pdf",
"png"))
{
par(ask = FALSE)
imageType <- match.arg(imageType)
indexes = start_ind:end_ind
temp2 = which(R[start_ind:end_ind, 4] == 0)
temp3 = indexes[temp2]
data_out_store = data_out
data_in_store = data_in
data_in = data_in[temp2]
data_out = data_out[temp2]
ldata_out = log(data_out)
n2 = length(ldata_out)
R[temp3, 5] = ldata_out
lag_phase_attempted = 2
MAXREG = 0
filename2 <- sprintf("%s_%s.%s", Name, code, imageType)
NUM = data_in[which(data_in <= FirstHours)]
if (n2 == 3) {
NE = 0
}
else if (length(NUM) <= 1) {
NE = 1
}
else if ((max(diff(NUM)) > MaxDiffInFirstHours) | (length(NUM) <
NeededInFirstHours)) {
NE = 1
}
else {
NE = 0
}
LP_temp = which(data_out != 0)
LLP = LP_temp[length(LP_temp)]
LastPos = data_out[LLP]
k = NULL
tlag = data_in[1]
if (n2 > 3) {
Y = log(data_out[1:(n2 - 1)])
X = data_in[1:(n2 - 1)]
cf = coef(lm(Y ~ X))
x_CI = try(confint(nls(Y ~ beta * (X - x0) - log(DT),
start = c(beta = cf[[2]], x0 = -cf[[1]]/cf[[2]]))))
if (class(x_CI) == "try-error") {
x_CI = array(NA, c(2, 2))
}
else if (is.na(x_CI[2, 1])) {
x_CI[2, 1] = NA
}
}
prof_type = -9999
if (n2 < 3) {
xx <- sprintf("TOBIT: There are not enough time points - no model fitted: data set %s",
code)
R[start_ind, 56] = 1
}
else if (data_out[1] < Threshold2) {
xx <- sprintf("TOBIT: Baseline level is too low - no model fitted: data set %s",
code)
R[start_ind, 56] = 2
}
else if (n2 == 3) {
ldata_out = log(data_out)
R[temp3, 5] = ldata_out
Linear = fitLinearAndStorage(data_in, ldata_out, start_ind,
temp3, R)
glm.linear = Linear$glm.linear
R = Linear$R
AIC1 = Linear$AIC1
mod = glm.linear
LinPart = seq(1, n2)
k = R[start_ind, 21]
alin = mod$coeff[[1]]
R[start_ind, 53] = R[start_ind, 22]
xx <- sprintf("TOBIT: Fitting a linear model since it has the lowest AIC n = 3: data set %s",
code)
prof_type = 1
R[start_ind, 44] = R[start_ind, 23]
}
else if ((data_in[n2] > x_CI[[2, 1]]) & (ldata_out[n2 - 1] >
Threshold3)) {
xx <- sprintf("TOBIT: zero is not informative- no model is fitted: data set %s",
code)
R[start_ind, 56] = 3.1
}
else if (NE == 1) {
ldata_out = log(data_out)
R[temp3, 5] = ldata_out
Linear = try(fitLinearAndStorage_tobit(data_in, ldata_out,
start_ind, temp3, R, DT, MaxIters))
glm.linear = Linear$glm.linear
R = Linear$R
AIC1 = Linear$AIC1
alin = Linear$a1
mod = glm.linear
LinPart = seq(1, n2)
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
k = R[start_ind, 21]
R[start_ind, 53] = R[start_ind, 22]
xx <- sprintf("TOBIT: Fitting a linear model since it has the lowest AIC n = 4: data set %s",
code)
prof_type = 1
R[start_ind, 44] = R[start_ind, 23]
}
else if (n2 == 4) {
ldata_out = log(data_out)
R[temp3, 5] = ldata_out
Linear = try(fitLinearAndStorage_tobit(data_in, ldata_out,
start_ind, temp3, R, DT, MaxIters))
glm.linear = Linear$glm.linear
R = Linear$R
AIC1 = Linear$AIC1
alin = Linear$a1
mod = glm.linear
LinPart = seq(1, n2)
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
k = R[start_ind, 21]
R[start_ind, 53] = R[start_ind, 22]
xx <- sprintf("TOBIT: Fitting a linear model since it has the lowest AIC n = 4: data set %s",
code)
prof_type = 1
R[start_ind, 44] = R[start_ind, 23]
}
else if (n2 == 5) {
Linear = try(fitLinearAndStorage_tobit(data_in, ldata_out,
start_ind, temp3, R, DT, MaxIters))
glm.linear = Linear$glm.linear
R = Linear$R
AIC1 = Linear$AIC1
alin = Linear$a1
Quad = try(fitQuadraticAndStorage_tobit(data_in, ldata_out,
start_ind, temp3, R, DT, MaxIters))
glm.quad = Quad$glm.quad
R = Quad$R
AIC2 = Quad$AIC2
scaledQ = Quad$scaled
ind = which(c(AIC1, AIC2) == min(c(AIC1, AIC2)))
ind = which(c(AIC1, AIC2) == min(c(AIC1, AIC2)))
if (data_out[2] > (1 + fact) * data_out[1]) {
data_in_mod = data_in[2:n2]
ldata_out_mod = ldata_out[2:n2]
tlag = data_in_mod[1]
m1 <- lm(ldata_out_mod ~ poly(data_in_mod, 1, raw = TRUE))
MaxReg <- try(tobit(ldata_out_mod ~ poly(data_in_mod,
1, raw = TRUE), left = log(DT), iter.max = MaxIters))
if (class(MaxReg)[1] == "try-error") {
yy <- sprintf("Fixing error in tobit fitting, data set: %s",
code)
print(yy)
options(warn = 1)
MaxReg <- tobit(ldata_out_mod ~ poly(data_in_mod,
1, raw = TRUE), left = log(DT), init = coef(m1),
iter.max = MaxIters)
options(warn = 2)
}
AICMaxReg = AIC(MaxReg)
ss_1 = sum(residuals(glm.linear)[2:(n2 - 1)]^2)
ss_2 = sum(residuals(glm.quad)[2:(n2 - 1)]^2)
ss_3 = sum(residuals(MaxReg)[1:3]^2)
ind = which(c(ss_1, ss_2, ss_3) == min(c(ss_1, ss_2,
ss_3)))
}
if (ind == 1) {
mod = glm.linear
LinPart = seq(1, n2)
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
tlag = data_in[1]
k = R[start_ind, 21]
R[start_ind, 53] = R[start_ind, 22]
xx <- sprintf("TOBIT: Fitting a linear model since it has the lowest SS n = 5: data set %s",
code)
prof_type = 1
R[start_ind, 44] = R[start_ind, 23]
}
else if (ind == 2) {
mod = glm.quad
mu = mean(data_in)
sd = sd(data_in)
a1 = coef(mod)[[1]]
b1 = coef(mod)[[2]]
c1 = coef(mod)[[3]]
if (scaledQ == 1) {
c = a1 - b1 * mu/sd + c1 * mu^2/sd^2
b = b1/sd - 2 * mu * c1/sd^2
a = c1/sd^2
}
else {
c = coef(mod)[[1]]
b = coef(mod)[[2]]
a = coef(mod)[[3]]
}
mod.fitted <- fitted.values(mod)
slopes = diff(mod.fitted)/diff(data_in)
c_points = c(data_in[1], data_in[length(data_in) -
1])
c_devs = 2 * a * c_points + b
ind_max = 2
fastest_rate = c_devs[ind_max]
ratios = fastest_rate/slopes
LinPart = c(which((ratios <= 5) & (ratios >= 0)),
n2)
ldattemp = mod.fitted[LinPart]
if (a < 0) {
xx <- sprintf("TOBIT: Fitting a quadratic (convex) model since it has the lowest SS and n = 5: data set %s",
code)
prof_type = 3
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
ldattemp = mod.fitted[LinPart]
tlag = data_in[LinPart[1]]
dattemp = data_in[LinPart]
mod <- glm(ldattemp ~ dattemp)
k = mod$coeff[[2]]
alin = mod$coeff[[1]]
}
else if (a >= 0) {
xx <- sprintf("TOBIT: Fitting a linear model since we have identified a concave quadratic (quadratic has the lowest SS) and n = 5: data set %s",
code)
prof_type = 2
LinPart = seq(1, n2)
ldattemp = ldata_out[LinPart]
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
tlag = data_in[LinPart[1]]
mod <- glm.linear
k = R[start_ind, 21]
}
if (tlag == data_in[1]) {
LinPart = seq(1, n2)
LinPart = LinPart[1:(length(LinPart) - 1)]
mod <- glm.linear
k = R[start_ind, 21]
to_include_R2 = LinPart
}
else {
to_include_R2 = seq(1, length(LinPart))
}
temp3 = indexes[temp2]
R[start_ind, 44] = 1 - sum((fitted(mod)[to_include_R2] -
ldata_out[LinPart])^2)/sum((ldata_out[LinPart] -
mean(ldata_out[LinPart]))^2)
SE_temp = R[start_ind, 22]
if (tlag > data_in[1]) {
SE_temp = sqrt(diag(vcov(mod)))[[2]]
}
R[start_ind, 53] = SE_temp
R[start_ind, 50] <- mod.fitted[which(data_in == tlag)[1]] -
mod.fitted[1]
R[start_ind, 50] <- mod.fitted[which(data_in == tlag)[1]] -
mod.fitted[1]
glm.quad.fitted <- fitted.values(glm.quad)
R[start_ind, 51] <- max(glm.quad.fitted) - glm.quad.fitted[which(data_in ==
tlag)[1]]
}
else if (ind == 3) {
mod <- MaxReg
AIC1 = AIC(mod)
LinPart = 1:(n2 - 1)
if (LinPart[length(LinPart)] == (n2 - 1)) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
k = coef(mod)[[2]]
alin = mod$coeff[[1]]
SEs2 = sqrt(diag(vcov(mod)))
R[start_ind, 53] = SEs2[[2]]
xx <- sprintf("TOBIT: Fitting a linear model from max value since we only have 5 data points and an initial increase in parasite level (also - lowest SS): data set %s",
code)
prof_type = 5
R[start_ind, 44] = 1 - sum(residuals(mod)^2)/sum((ldata_out_mod -
mean(ldata_out_mod))^2)
MAXREG = 1
}
}
else {
Linear = fitLinearAndStorage_tobit(data_in, ldata_out,
start_ind, temp3, R, DT, MaxIters)
glm.linear = Linear$glm.linear
R = Linear$R
AIC1 = Linear$AIC1
alin = Linear$a1
Quad = fitQuadraticAndStorage_tobit(data_in, ldata_out,
start_ind, temp3, R, DT, MaxIters)
glm.quad = Quad$glm.quad
R = Quad$R
AIC2 = Quad$AIC2
scaledQ = Quad$scaled
Cubic = fitCubicAndStorage_tobit(data_in, ldata_out,
start_ind, temp3, R, DT, MaxIters)
glm.cubic = Cubic$glm.cubic
R = Cubic$R
AIC3 = Cubic$AIC3
scaledC = Cubic$scaled
ind = which(c(AIC1, AIC2, AIC3) == min(c(AIC1, AIC2,
AIC3)))
if (ind == 1) {
mod = glm.linear
LinPart = seq(1, n2)
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
k = R[start_ind, 21]
R[start_ind, 53] = R[start_ind, 22]
xx <- sprintf("TOBIT: Fitting a linear model since it has the lowest AIC: data set %s",
code)
prof_type = 1
R[start_ind, 44] = R[start_ind, 23]
}
else if (ind == 2) {
mod = glm.quad
mu = mean(data_in)
sd = sd(data_in)
a1 = coef(mod)[[1]]
b1 = coef(mod)[[2]]
c1 = coef(mod)[[3]]
c = a1 - b1 * mu/sd + c1 * mu^2/sd^2
b = b1/sd - 2 * mu * c1/sd^2
a = c1/sd^2
if (scaledQ == 1) {
c = a1 - b1 * mu/sd + c1 * mu^2/sd^2
b = b1/sd - 2 * mu * c1/sd^2
a = c1/sd^2
}
else {
c = coef(mod)[[1]]
b = coef(mod)[[2]]
a = coef(mod)[[3]]
}
mod.fitted <- fitted.values(mod)
slopes = diff(mod.fitted)/diff(data_in)
c_points = c(data_in[1], data_in[length(data_in) -
1])
c_devs = 2 * a * c_points + b
ind_max = 2
fastest_rate = c_devs[ind_max]
ratios = fastest_rate/slopes
LinPart = c(which((ratios <= 5) & (ratios >= 0)),
n2)
ldattemp = mod.fitted[LinPart]
if (a < 0) {
xx <- sprintf("TOBIT: Fitting a quadratic (convex) model since it has the lowest AIC: data set %s",
code)
prof_type = 3
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
ldattemp = mod.fitted[LinPart]
tlag = data_in[LinPart[1]]
dattemp = data_in[LinPart]
mod <- glm(ldattemp ~ dattemp)
k = mod$coeff[[2]]
alin = mod$coeff[[1]]
}
else if (a >= 0) {
xx <- sprintf("TOBIT: Fitting a linear model since we have identified a concave quadratic (quadratic has the lowest AIC): data set %s",
code)
prof_type = 2
LinPart = seq(1, n2)
ldattemp = ldata_out[LinPart]
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
tlag = data_in[LinPart[1]]
mod <- glm.linear
k = R[start_ind, 21]
}
if (tlag == data_in[1]) {
LinPart = seq(1, n2)
LinPart = LinPart[1:(length(LinPart) - 1)]
mod <- glm.linear
k = R[start_ind, 21]
to_include_R2 = LinPart
}
else {
to_include_R2 = seq(1, length(LinPart))
}
temp3 = indexes[temp2]
R[start_ind, 44] = 1 - sum((fitted(mod)[to_include_R2] -
ldata_out[LinPart])^2)/sum((ldata_out[LinPart] -
mean(ldata_out[LinPart]))^2)
SE_temp = R[start_ind, 22]
if (tlag > data_in[1]) {
SE_temp = sqrt(diag(vcov(mod)))[[2]]
}
R[start_ind, 53] = SE_temp
R[start_ind, 50] <- mod.fitted[which(data_in == tlag)[1]] -
mod.fitted[1]
R[start_ind, 50] <- mod.fitted[which(data_in == tlag)[1]] -
mod.fitted[1]
glm.quad.fitted <- fitted.values(glm.quad)
R[start_ind, 51] <- max(glm.quad.fitted) - glm.quad.fitted[which(data_in ==
tlag)[1]]
}
else if (ind == 3) {
mod = glm.cubic
a1 = coef(mod)[[1]]
b1 = coef(mod)[[2]]
c1 = coef(mod)[[3]]
d1 = coef(mod)[[4]]
mu = mean(data_in)
sd = sd(data_in)
if (scaledC == 1) {
d = a1 - b1 * mu/sd + c1 * mu^2/sd^2 - d1 * mu^3/sd^3
c = b1/sd - 2 * mu * c1/sd^2 + 3 * mu^2 * d1/sd^3
b = c1/sd^2 - 3 * mu * d1/sd^3
a = d1/sd^3
}
else {
d = coef(mod)[[1]]
c = coef(mod)[[2]]
b = coef(mod)[[3]]
a = coef(mod)[[4]]
}
PoI = -b/(3 * a)
c_points = c(data_in[1], PoI, data_in[length(data_in) -
1])
c_devs = 3 * a * c_points^2 + 2 * b * c_points +
c
ind_min = which(abs(c_devs) == min(abs(c_devs)))
slowest_rate = c_devs[ind_min]
ind_max = which(abs(c_devs) == max(abs(c_devs)))
fastest_rate = c_devs[which(c_devs == min(c_devs))]
ind_max_neg = which(c_devs == min(c_devs))
fastest_rate = c_devs[ind_max_neg]
mod.fitted <- fitted.values(mod)
slopes = diff(mod.fitted)/diff(data_in)
ratios = fastest_rate/slopes
FastParts = which((ratios <= 5) & (ratios >= 0))
if (((PoI < data_in[1]) | (PoI > data_in[length(data_in) -
1])) & (b < 0)) {
xx <- sprintf("TOBIT: Fitting a cubic model for cubic case (d) where b <0: data set %s",
code)
fastest_rate = min(c(c_devs[1], c_devs[3]))
ratios = fastest_rate/slopes
prof_type = 7
LinPart = which((ratios <= 5) & (ratios >= 0))
temp = seq(1, n2)
LinPart = c(LinPart, temp[LinPart[length(LinPart)] +
1])
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
ldattemp = mod.fitted[LinPart]
tlag = data_in[LinPart[1]]
dattemp = data_in[LinPart]
mod <- glm(ldattemp ~ dattemp)
k = mod$coeff[[2]]
alin = mod$coeff[[1]]
}
else if (((PoI < data_in[1]) | (PoI > data_in[length(data_in) -
1])) & (b > 0)) {
xx <- sprintf("TOBIT: Fitting a linear model for cubic case (d) where b >0: data set %s",
code)
prof_type = 6
LinPart = seq(1, n2)
ldattemp = ldata_out[LinPart]
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
tlag = data_in[LinPart[1]]
mod <- glm.linear
k = R[start_ind, 21]
}
else if (ind_max_neg == 2) {
xx <- sprintf("TOBIT: Fitting a cubic model for cubic case (a): data set %s",
code)
prof_type = 4
LinPart = which((ratios <= 5) & (ratios >= 0) &
(slopes < 0))
temp = seq(1, n2)
LinPart = c(LinPart, temp[LinPart[length(LinPart)] +
1])
ldattemp = mod.fitted[LinPart]
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
ldattemp = mod.fitted[LinPart]
tlag = data_in[LinPart[1]]
dattemp = data_in[LinPart]
mod <- glm(ldattemp ~ dattemp)
k = mod$coeff[[2]]
alin = mod$coeff[[1]]
}
else if ((ind_max_neg != 2) & (FastParts[1] == 1)) {
xx <- sprintf("TOBIT: Fitting a linear model for cubic case (biii) where there is no lag phase: data set %s",
code)
prof_type = 9
LinPart = seq(1, n2)
ldattemp = ldata_out[LinPart]
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
tlag = data_in[LinPart[1]]
mod <- glm.linear
k = R[start_ind, 21]
}
else if ((ind_max_neg != 2) & (FastParts[1] != 1)) {
xx <- sprintf("TOBIT: Fitting a linear model for cubic case (bii) where there is a lag phase: data set %s",
code)
prof_type = 8
LinPart = which((ratios <= 5) & (ratios >= 0) &
(data_in[2:n2] > PoI))
temp = seq(1, n2)
LinPart = c(LinPart, temp[LinPart[length(LinPart)] +
1])
ldattemp = mod.fitted[LinPart]
if (LinPart[length(LinPart)] == n2) {
LinPart = LinPart[1:(length(LinPart) - 1)]
}
ldattemp = mod.fitted[LinPart]
tlag = data_in[LinPart[1]]
dattemp = data_in[LinPart]
mod <- glm(ldattemp ~ dattemp)
k = mod$coeff[[2]]
alin = mod$coeff[[1]]
}
if (tlag == data_in[1]) {
LinPart = seq(1, n2)
LinPart = LinPart[1:(length(LinPart) - 1)]
mod <- glm.linear
k = R[start_ind, 21]
to_include_R2 = LinPart
}
else {
to_include_R2 = seq(1, length(LinPart))
}
temp3 = indexes[temp2]
R[start_ind, 44] = 1 - sum((fitted(mod)[to_include_R2] -
ldata_out[LinPart])^2)/sum((ldata_out[LinPart] -
mean(ldata_out[LinPart]))^2)
SE_temp = R[start_ind, 22]
if (tlag > data_in[1]) {
SE_temp = sqrt(diag(vcov(mod)))[[2]]
}
R[start_ind, 53] = SE_temp
R[start_ind, 50] <- mod.fitted[which(data_in == tlag)[1]] -
mod.fitted[1]
glm.cubic.fitted <- fitted.values(glm.cubic)
R[start_ind, 51] <- max(glm.cubic.fitted) - glm.cubic.fitted[which(data_in ==
tlag)[1]]
R[start_ind, 51] <- max(glm.cubic.fitted) - glm.cubic.fitted[which(data_in ==
tlag)[1]]
}
}
if ((plotting == TRUE) & all(ldata_out != -Inf)) {
if (imageType == "png") {
Cairo(file = filename2, type = "png")
}
else if (imageType == "pdf") {
pdf(file = filename2)
}
filename <- sprintf("patient id: %s", code)
par(mar = c(5.1, 5.1, 4.1, 2.1))
x0 = data_in[1]
xend = data_in[length(data_in)] * 1.05
y0 = 0
yend = max(ldata_out) * 1.05
if (yend == -Inf) {
yend = 10
}
xlim = c(x0, xend)
ylim = c(y0, yend)
plot(data_in, ldata_out, main = filename, xlab = "Time (hours)",
ylab = expression(log[e](parasitaemia)), cex = 2,
cex.main = 2, cex.lab = 2, cex.axis = 2, xlim = xlim,
ylim = ylim, type = "p", pch = 19)
if (R[start_ind, 56] == 0) {
if (length(fitted(mod)) == length(LinPart)) {
lines(data_in[LinPart], fitted(mod), lty = 1,
col = "blue", lwd = 3)
}
else if (prof_type == 5) {
lines(data_in[LinPart + 1], fitted(mod)[LinPart],
lty = 1, col = "blue", lwd = 3)
}
else {
lines(data_in[LinPart], fitted(mod)[LinPart],
lty = 1, col = "blue", lwd = 3)
}
}
out_temp = floor(as.numeric(R[start_ind:end_ind, 4]))
outliers = which(out_temp == 2)
if (length(outliers) > 0) {
points(data_in_store[outliers], log(data_out_store[outliers]),
col = "red", pch = 19, cex = 2)
}
if ((tlag > 0) & (R[start_ind, 56] == 0)) {
t_temp = which(data_in < tlag)
points(data_in[t_temp], log(data_out[t_temp]), col = "gray",
pch = 19, cex = 2)
}
points(data_in_store[ind_DL], log(data_out_store[ind_DL]),
col = "green", pch = 19, cex = 2)
dev.off()
}
print(xx)
if (is.null(k) == FALSE) {
mod_lin <- lm(ldata_out[-ind_DL] ~ data_in[-ind_DL])
R[start_ind, 70] = 1 - sum(residuals(mod_lin)^2)/sum((ldata_out[-ind_DL] -
mean(ldata_out[-ind_DL]))^2)
if (MAXREG == 1) {
MSRE = sqrt(mean((fitted(mod)[LinPart] - ldata_out[LinPart +
1])^2))
}
else if (tlag == data_in[1]) {
MSRE = (mean(residuals(mod)^2))^(1/2)
}
else {
MSRE = sqrt(mean((fitted(mod) - ldata_out[LinPart])^2))
}
R[start_ind, 55] = MSRE
R[temp3[LinPart], 52] = 1
if (MAXREG == 1) {
R[temp3[LinPart + 1], 45] <- fitted.values(mod)[LinPart]
R[temp3[LinPart + 1], 46] <- residuals(mod)[LinPart]
}
else if (length(fitted(mod)) == length(LinPart)) {
R[temp3[LinPart], 45] <- fitted.values(mod)
R[temp3[LinPart], 46] <- residuals(mod)
}
else {
R[temp3[LinPart], 45] <- fitted.values(mod)[LinPart]
R[temp3[LinPart], 46] <- residuals(mod)[LinPart]
}
R[start_ind, 49] <- tlag - data_in[1]
R[start_ind, 42] <- prof_type
R[start_ind, 43] <- k
knum = as.numeric(k)
int = alin + knum * as.numeric(R[start_ind, 49])
R[start_ind, 67] = int
PC50 = (log(data_out[1] * (1 - 0.5)) - alin)/knum
PC90 = (log(data_out[1] * (1 - 0.9)) - alin)/knum
PC95 = (log(data_out[1] * (1 - 0.95)) - alin)/knum
PC99 = (log(data_out[1] * (1 - 0.99)) - alin)/knum
if (PC50 > 0) {
R[start_ind, 59] = PC50
}
if (PC90 > 0) {
R[start_ind, 60] = PC90
}
if (PC95 > 0) {
R[start_ind, 61] = PC95
}
if (PC99 > 0) {
R[start_ind, 62] = PC99
}
lag_phase_attempted = 1
if ((NE == 0) | (n2 == 3)) {
lag_phase_attempted = 0
}
R[start_ind, 50] <- data_out[which(data_in == tlag)[1]] -
max(data_out[which(data_in <= tlag)])
R[start_ind, 51] <- data_out[which(data_in == tlag)[1]] -
data_out[1]
if (k > 0) {
xx <- sprintf("TOBIT: NOTE: Model is fitted, but parasites not cleared: data set %s",
code)
print(xx)
}
}
R[start_ind, 57] <- lag_phase_attempted
return(list(R = R))
}
<bytecode: 0x9888b28>
<environment: namespace:bhrcr>
--- function search by body ---
Function lagReg_tobit in namespace bhrcr has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(x_CI) == "try-error") { : the condition has length > 1
Calls: calculatePCE
Execution halted
Flavor: r-devel-linux-x86_64-fedora-gcc