It appears to be an issue with argument handling in the logistic_model
function. Unpacking the pieces seems to run as desired.
# as in OP
library(tidyverse)
library(lubridate)
#>
#> Attaching package: 'lubridate'
#> The following objects are masked from 'package:dplyr':
#>
#> intersect, setdiff, union
#> The following objects are masked from 'package:base':
#>
#> date, intersect, setdiff, union
library(car)
#> Loading required package: carData
#>
#> Attaching package: 'car'
#> The following object is masked from 'package:dplyr':
#>
#> recode
#> The following object is masked from 'package:purrr':
#>
#> some
indata <- tribble(
~Date, ~Cases, ~Days,
"2020-03-10", 21, 0,
"2020-03-11", 23, 1,
"2020-03-12", 38, 2,
"2020-03-13", 51, 3,
"2020-03-14", 56, 4,
"2020-03-15", 57, 5,
"2020-03-16", 64, 6,
"2020-03-17", 83, 7,
"2020-03-18", 143, 8,
"2020-03-19", 194, 9,
"2020-03-20", 304, 10,
"2020-03-21", 334, 11,
"2020-03-22", 352, 12,
"2020-03-23", 410, 13,
"2020-03-24", 974, 14,
"2020-03-25", 1396, 15,
"2020-03-26", 1731, 16,
"2020-03-27", 2052, 17,
"2020-03-28", 2371, 18,
"2020-03-29", 2877, 19,
"2020-03-30", 3266, 20,
"2020-03-31", 3997, 21,
"2020-04-01", 4669, 22,
"2020-04-02", 5330, 23,
"2020-04-03", 6110, 24,
"2020-04-04", 6812, 25,
"2020-04-05", 7276, 26
)
indata$Date <- as_date(indata$Date)
# Unpack function to inspect steps
# logistic_model <- function(df.,
# indep="Cases", # independent variable
# r=0.24,
# projection=10){
df. <- indata # avoid potential namespace confict with df.()
df.
#> # A tibble: 27 x 3
#> Date Cases Days
#> <date> <dbl> <dbl>
#> 1 2020-03-10 21 0
#> 2 2020-03-11 23 1
#> 3 2020-03-12 38 2
#> 4 2020-03-13 51 3
#> 5 2020-03-14 56 4
#> 6 2020-03-15 57 5
#> 7 2020-03-16 64 6
#> 8 2020-03-17 83 7
#> 9 2020-03-18 143 8
#> 10 2020-03-19 194 9
#> # … with 17 more rows
r <- 0.24
projection <- 10
indep <- df.$Cases
Asym <- max(indep)
Asym
#> [1] 7276
Asym <- max(indep)*5
Asym
#> [1] 36380
xmid <- max(indep)*2
xmid
#> [1] 14552
scal <- 1/r
scal
#> [1] 4.166667
## using a selfStart model
my_model <- nls(Cases ~ SSlogis(Days, Asym, xmid, scal), data=df.)
my_model
#> Nonlinear regression model
#> model: Cases ~ SSlogis(Days, Asym, xmid, scal)
#> data: df.
#> Asym xmid scal
#> 10130.468 22.497 3.706
#> residual sum-of-squares: 357740
#>
#> Number of iterations to convergence: 0
#> Achieved convergence tolerance: 1.495e-06
coeffs <- coef(my_model)
coeffs
#> Asym xmid scal
#> 10130.467944 22.496894 3.705809
dayseq <- df.$Days
dayseq
#> [1] 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
#> [26] 25 26
dayseq <- c(dayseq,(dayseq[length(dayseq)]+1):
(dayseq[length(dayseq)]-1+projection))
dayseq
#> [1] 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
#> [26] 25 26 27 28 29 30 31 32 33 34 35
dateseq <- df.$Date
dateseq
#> [1] "2020-03-10" "2020-03-11" "2020-03-12" "2020-03-13" "2020-03-14"
#> [6] "2020-03-15" "2020-03-16" "2020-03-17" "2020-03-18" "2020-03-19"
#> [11] "2020-03-20" "2020-03-21" "2020-03-22" "2020-03-23" "2020-03-24"
#> [16] "2020-03-25" "2020-03-26" "2020-03-27" "2020-03-28" "2020-03-29"
#> [21] "2020-03-30" "2020-03-31" "2020-04-01" "2020-04-02" "2020-04-03"
#> [26] "2020-04-04" "2020-04-05"
dateseq <- as_date(c(dateseq,(dateseq[length(dateseq)]+1):
(dateseq[length(dateseq)]-1+projection)))
dateseq
#> [1] "2020-03-10" "2020-03-11" "2020-03-12" "2020-03-13" "2020-03-14"
#> [6] "2020-03-15" "2020-03-16" "2020-03-17" "2020-03-18" "2020-03-19"
#> [11] "2020-03-20" "2020-03-21" "2020-03-22" "2020-03-23" "2020-03-24"
#> [16] "2020-03-25" "2020-03-26" "2020-03-27" "2020-03-28" "2020-03-29"
#> [21] "2020-03-30" "2020-03-31" "2020-04-01" "2020-04-02" "2020-04-03"
#> [26] "2020-04-04" "2020-04-05" "2020-04-06" "2020-04-07" "2020-04-08"
#> [31] "2020-04-09" "2020-04-10" "2020-04-11" "2020-04-12" "2020-04-13"
#> [36] "2020-04-14"
foo <- tibble(Days=dayseq)
foo
#> # A tibble: 36 x 1
#> Days
#> <dbl>
#> 1 0
#> 2 1
#> 3 2
#> 4 3
#> 5 4
#> 6 5
#> 7 6
#> 8 7
#> 9 8
#> 10 9
#> # … with 26 more rows
predict2 <- function(x) {stats::predict(x, foo)}
print(summary(foo))
#> Days
#> Min. : 0.00
#> 1st Qu.: 8.75
#> Median :17.50
#> Mean :17.50
#> 3rd Qu.:26.25
#> Max. :35.00
f.boot <- car::Boot(my_model,f=predict2)
#> Loading required namespace: boot
f.boot
#>
#> ORDINARY NONPARAMETRIC BOOTSTRAP
#>
#>
#> Call:
#> boot::boot(data = data.frame(update(object, model = TRUE)$model),
#> statistic = boot.f, R = R, .fn = f, parallel = p_type, ncpus = ncores,
#> cl = cl2)
#>
#>
#> Bootstrap Statistics :
#> original bias std. error
#> t1* 23.34271 -0.4866764 6.170080
#> t2* 30.55162 -0.6494386 7.634477
#> t3* 39.97806 -0.8572582 9.413758
#> t4* 52.29787 -1.1194215 11.561217
#> t5* 68.38849 -1.4456372 14.132285
#> t6* 89.38583 -1.8450955 17.180618
#> t7* 116.75519 -2.3248455 20.751736
#> t8* 152.37769 -2.8872076 24.873422
#> t9* 198.65316 -3.5259009 29.542128
#> t10* 258.61775 -4.2206262 34.705157
#> t11* 336.07050 -4.9301107 40.239620
#> t12* 435.69562 -5.5842901 45.931675
#> t13* 563.15564 -6.0776290 51.463497
#> t14* 725.11438 -6.2677814 56.419930
#> t15* 929.12766 -5.9866532 60.328484
#> t16* 1183.31907 -5.0732167 62.738629
#> t17* 1495.74780 -3.4361798 63.324491
#> t18* 1873.39421 -1.1456286 61.971226
#> t19* 2320.76197 1.4664003 58.829197
#> t20* 2838.23972 3.7410084 54.471486
#> t21* 3420.56561 4.7221588 50.586974
#> t22* 4055.91561 3.3230391 51.475973
#> t23* 4726.15538 -1.3874649 63.485987
#> t24* 5408.53856 -9.8738393 89.480732
#> t25* 6078.62835 -21.9394712 127.932482
#> t26* 6713.69899 -36.7506952 176.224488
#> t27* 7295.65160 -53.0699325 231.491267
#> t28* 7812.70357 -69.5831710 290.621872
#> t29* 8259.63192 -85.1818215 350.520681
#> t30* 8636.85592 -99.1149886 408.479233
#> t31* 8948.89958 -111.0078632 462.444355
#> t32* 9202.75391 -120.7925148 511.108201
#> t33* 9406.48122 -128.6062318 553.848349
#> t34* 9568.20315 -134.6962370 590.584056
#> t35* 9695.47071 -139.3483623 621.609453
#> t36* 9794.94166 -142.8420482 647.442477
Cases <- predict(my_model, data.frame(Days=dayseq))
Cases
#> [1] 23.34271 30.55162 39.97806 52.29787 68.38849 89.38583
#> [7] 116.75519 152.37769 198.65316 258.61775 336.07050 435.69562
#> [13] 563.15564 725.11438 929.12766 1183.31907 1495.74780 1873.39421
#> [19] 2320.76197 2838.23972 3420.56561 4055.91561 4726.15538 5408.53856
#> [25] 6078.62835 6713.69899 7295.65160 7812.70357 8259.63192 8636.85592
#> [31] 8948.89958 9202.75391 9406.48122 9568.20315 9695.47071 9794.94166
#> attr(,"gradient")
#> Asym xmid scal
#> [1,] 0.002304208 -6.284438 38.15101
#> [2,] 0.003015816 -8.219390 47.67957
#> [3,] 0.003946319 -10.745372 59.43284
#> [4,] 0.005162434 -14.039551 73.86448
#> [5,] 0.006750773 -18.329821 91.49009
#> [6,] 0.008823465 -23.907639 112.87939
#> [7,] 0.011525152 -31.142883 138.63664
#> [8,] 0.015041525 -40.500121 169.36279
#> [9,] 0.019609475 -52.554703 205.59075
#> [10,] 0.025528707 -68.005554 247.68245
#> [11,] 0.033174233 -87.679007 295.67508
#> [12,] 0.043008440 -112.514452 349.06464
#> [13,] 0.055590289 -143.517840 406.52168
#> [14,] 0.071577580 -181.664114 465.55152
#> [15,] 0.091716164 -227.726715 522.14509
#> [16,] 0.116807938 -282.016183 570.52203
#> [17,] 0.147648441 -344.028278 603.13834
#> [18,] 0.184926720 -412.043283 611.19138
#> [19,] 0.229087342 -482.783908 585.84463
#> [20,] 0.280168669 -551.311217 520.23111
#> [21,] 0.337651294 -611.366485 411.92558
#> [22,] 0.400368041 -656.282299 265.09332
#> [23,] 0.466528833 -680.355590 91.22562
#> [24,] 0.533888325 -680.278775 -92.35562
#> [25,] 0.600034311 -656.062687 -266.10432
#> [26,] 0.662723482 -611.033443 -412.72542
#> [27,] 0.720169260 -550.904750 -520.77100
#> [28,] 0.771208558 -482.345390 -586.12103
#> [29,] 0.815325803 -411.608108 -611.23583
#> [30,] 0.852562386 -343.622014 -603.00209
#> [31,] 0.883364878 -281.654045 -570.26153
#> [32,] 0.908423378 -227.415172 -521.81197
#> [33,] 0.928533733 -181.403348 -465.18733
#> [34,] 0.944497648 -143.304155 -406.15662
#> [35,] 0.957060499 -112.342194 -348.71853
#> [36,] 0.966879488 -87.541887 -295.35942
confidence <- confint(f.boot)
confidence
#> Bootstrap bca confidence intervals
#>
#> 2.5 % 97.5 %
#> V1 12.29421 36.78035
#> V2 16.62235 46.77471
#> V3 22.50012 59.62103
#> V4 30.35387 76.08839
#> V5 41.50583 97.16222
#> V6 56.45748 123.92827
#> V7 76.58277 157.73101
#> V8 103.19305 200.37641
#> V9 138.81629 254.28667
#> V10 186.76832 322.77210
#> V11 252.09683 409.49938
#> V12 339.23453 518.00073
#> V13 453.94899 655.36019
#> V14 601.14273 825.64160
#> V15 793.03875 1036.14563
#> V16 1038.62081 1288.26522
#> V17 1346.21539 1606.07691
#> V18 1725.10013 1976.35854
#> V19 2180.92446 2421.26546
#> V20 2716.17593 2935.47340
#> V21 3308.69274 3515.97235
#> V22 3974.49671 4167.46224
#> V23 4650.81530 4848.08154
#> V24 5297.29771 5530.97005
#> V25 5867.97492 6179.78124
#> V26 6421.29890 6796.79802
#> V27 7071.33555 7490.84967
#> V28 7477.23413 8129.73858
#> V29 7754.22767 8669.90743
#> V30 7955.41737 9126.02168
#> V31 8151.33134 9552.85961
#> V32 8284.31206 9913.10172
#> V33 8360.49797 10191.01474
#> V34 8422.28792 10417.73427
#> V35 8501.58178 10634.02729
#> V36 8551.54545 10792.10457
preds <- tibble(dayseq, dateseq,Cases, confidence["2.5 %"], confidence["97.5 %"])
preds
#> # A tibble: 36 x 5
#> dayseq dateseq Cases `confidence["2.5 %"]` `confidence["97.5 %"]`
#> <dbl> <date> <dbl> <dbl> <dbl>
#> 1 0 2020-03-10 23.3 NA NA
#> 2 1 2020-03-11 30.6 NA NA
#> 3 2 2020-03-12 40.0 NA NA
#> 4 3 2020-03-13 52.3 NA NA
#> 5 4 2020-03-14 68.4 NA NA
#> 6 5 2020-03-15 89.4 NA NA
#> 7 6 2020-03-16 117. NA NA
#> 8 7 2020-03-17 152. NA NA
#> 9 8 2020-03-18 199. NA NA
#> 10 9 2020-03-19 259. NA NA
#> # … with 26 more rows
names(preds) <- c("Days", "Date","Cases","SD_upper","SD_lower")
preds
#> # A tibble: 36 x 5
#> Days Date Cases SD_upper SD_lower
#> <dbl> <date> <dbl> <dbl> <dbl>
#> 1 0 2020-03-10 23.3 NA NA
#> 2 1 2020-03-11 30.6 NA NA
#> 3 2 2020-03-12 40.0 NA NA
#> 4 3 2020-03-13 52.3 NA NA
#> 5 4 2020-03-14 68.4 NA NA
#> 6 5 2020-03-15 89.4 NA NA
#> 7 6 2020-03-16 117. NA NA
#> 8 7 2020-03-17 152. NA NA
#> 9 8 2020-03-18 199. NA NA
#> 10 9 2020-03-19 259. NA NA
#> # … with 26 more rows
# return a tibble
tribble(~Line, ~r, ~K, ~xmid,
preds, 1/coeffs[["scal"]], coeffs[["Asym"]], coeffs[["xmid"]])
#> # A tibble: 1 x 4
#> Line r K xmid
#> <list> <dbl> <dbl> <dbl>
#> 1 <tibble [36 × 5]> 0.270 10130. 22.5
#}
Created on 2020-04-06 by the reprex package (v0.3.0)