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: 0x558118f7ac60>
<environment: namespace:bhrcr>
 --- function search by body ---
Function lagReg_tobit in namespace bhrcr has this body.
 ----------- END OF FAILURE REPORT -------------- 
Fatal error: the condition has length > 1