testing_bias_distribution/float.data.sim.R at main · bbolker/testing_bias_distribution · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
# Sys.setenv(LANG = "en")

library(shellpipes)
rpcall("const.data.Rout float.data.sim.R const.params.rda")
rpcall("float.data.Rout float.data.sim.R float.params.rda")
suppressPackageStartupMessages(library(dplyr))
library(tidyr)

suppressPackageStartupMessages(library(macpan2))
options(macpan2_verbose = FALSE)

#library(ggplot2); theme_set(theme_bw())
loadEnvironments()

#source("mle2_tidy.R")
#source("spec_trans_par.R")
#source("float.params.R")

# Set Seeds
set.seed(seed)

gamma <- deltat/D
t <- c(tmin:tmax)
pts <- length(t) ## number of time points

logit_trans <- function(x){
  log(x)-log(1-x)
}
logit_backtrans <- function(x){
  (1/(1 + exp(-x)))
}

## BMB: use self-naming list from tibble pkg
true_param <- tibble::lst(  log_h=log(h)
                          , log_I0=log(I0)
                          , log_W0=log(w0)
                          , log_WI=log(wI)
                          , log_alpha=log(alpha)
                          )
print(true_param)

tp_list <-tibble::lst(  beta
                      , gamma
                      , N
                      , h
                      , w0
                      , wI
                      , alpha
                      , I = I0
                      , R = 0
)

### import SIR from macpan-model-lib
mc_sir <- mp_tmb_library("starter_models","sir", package = "macpan2")

(mc_sir
  |> mp_tmb_update(
    default = tp_list
    )
  |> mp_tmb_insert_backtrans(  variables = c( "beta"
                                             ,"gamma"
                                             ,"I","h"
                                             ,"w0"
                                             ,"wI"
                                             ,"alpha")
                             , mp_log)
  |> mp_tmb_insert(
      phase = "during"
    , at = Inf
    , expressions = list(
        pY ~ I/N                          ## Prevalence based on SIR
      , pSus ~ S/N                        ## Susceptible proportion
      , b ~ w0+wI*pY*exp(-alpha*(1-pSus)) ## Concern floating baseline hazard
      , T_B ~ 1-exp(-b)
      , T_Y ~ 1-exp(-b-h)
      , T_prop ~ (1-pY)*T_B+pY*T_Y        ## Expected test proportion
      , T_pos ~ pY*T_Y/T_prop             ## Expected test positivity
      , OT ~ rbinom(N,T_prop)
      , OP ~ rbinom(OT,T_pos)
    )
  )
  ## |> mp_tmb_delete(phase = "before", at = Inf, default = c("beta","gamma","I","T_B","T_Y"))
)->sir

# sir |> mp_expand()

(sir
  |> mp_simulator(
      time_steps = tmax
    , outputs = c("pY","T_prop","T_pos","OT","OP","I","S","T_B")
  )
  |> mp_trajectory()
  |> dplyr::select(-c(row, col))

) -> dat_all

dat<- (dat_all|> filter(time>=tmin)
)


(dat
  |> pivot_wider(names_from = matrix,values_from = value)
) -> dat_fit

(dat_all
  |> pivot_wider(names_from = matrix,values_from = value)
) -> dat_pall

saveEnvironment()