R version 3.6.1 (2019-07-05) -- "Action of the Toes" Copyright (C) 2019 The R Foundation for Statistical Computing Platform: x86_64-pc-linux-gnu (64-bit) R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. Natural language support but running in an English locale R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > params <- + list(prefix = "decomp") > > ## ----setup,message=FALSE,cache=FALSE,purl=TRUE--------------------------- > library(plyr) > library(tidyverse) ── Attaching packages ─────────────────────────────────────── tidyverse 1.2.1 ── ✔ ggplot2 3.2.1 ✔ purrr 0.3.2 ✔ tibble 2.1.3 ✔ dplyr 0.8.3 ✔ tidyr 1.0.0 ✔ stringr 1.4.0 ✔ readr 1.3.1 ✔ forcats 0.4.0 ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ── ✖ dplyr::arrange() masks plyr::arrange() ✖ purrr::compact() masks plyr::compact() ✖ dplyr::count() masks plyr::count() ✖ dplyr::failwith() masks plyr::failwith() ✖ dplyr::filter() masks stats::filter() ✖ dplyr::id() masks plyr::id() ✖ dplyr::lag() masks stats::lag() ✖ dplyr::mutate() masks plyr::mutate() ✖ dplyr::rename() masks plyr::rename() ✖ dplyr::summarise() masks plyr::summarise() ✖ dplyr::summarize() masks plyr::summarize() > library(magrittr) Attaching package: ‘magrittr’ The following object is masked from ‘package:purrr’: set_names The following object is masked from ‘package:tidyr’: extract > library(tibble) > library(stringi) > library(pomp) Welcome to pomp version 2! For information on upgrading your pomp version < 2 code, see the ‘pomp version 2 upgrade guide’ at https://kingaa.github.io/pomp/. Attaching package: ‘pomp’ The following object is masked from ‘package:purrr’: map > library(panelPomp) > library(foreach) Attaching package: ‘foreach’ The following objects are masked from ‘package:purrr’: accumulate, when > library(iterators) > library(doRNG) Loading required package: rngtools Loading required package: pkgmaker Loading required package: registry Attaching package: ‘pkgmaker’ The following object is masked from ‘package:base’: isFALSE > library(aakmisc) ## available at https://kingaa.github.io/ setting ‘stringsAsFactors=FALSE’ > stopifnot(packageVersion("pomp")>="2.2") > stopifnot(packageVersion("panelPomp")>="0.9.1") > stopifnot(packageVersion("aakmisc")>="0.26.2") > options( + stringsAsFactors=FALSE, + keep.source=TRUE, + encoding="UTF-8" + ) > set.seed(407958184) > > > > > ## ----import_data,purl=TRUE----------------------------------------------- > read_csv("data.csv", + col_types="iiinnnn" + ) %>% + mutate( + mouseid=sprintf("%02d-%02d",box,mouse), + box=sprintf("%02d",box), + mouse=sprintf("%02d",mouse), + paba=as.factor(box), + rbc_density=rbc_density/1000 + ) %>% + mutate( + paba=recode( + paba, + "01"="0.05","02"="0.005","03"="0.0005","04"="0","05"="control" + ) + ) %>% + select( + day, + Pd=ama_density, + RBC=rbc_density, + Ter119=ter119_density, + CD71=cd71_density, + mouseid, + paba,box,mouse + ) %>% + arrange(mouseid,day) %>% + mutate( + paba=as.character(paba), + Ter119=ifelse(Ter119==0,NA,Ter119), + CD71=ifelse(CD71==0,NA,CD71) + ) %>% + mutate( + Eryth=(1-CD71/Ter119)*RBC, + Retic=CD71/Ter119*RBC + ) -> flow > > > > > > > ## ----beta-dose,purl=TRUE------------------------------------------------- > flow %>% + filter(day<=4) %>% + lm(log(Pd)~day+mouseid-1,data=.) -> fit1 > summary(fit1) Call: lm(formula = log(Pd) ~ day + mouseid - 1, data = .) Residuals: Min 1Q Median 3Q Max -0.63666 -0.28417 -0.02179 0.33097 0.76375 Coefficients: Estimate Std. Error t value Pr(>|t|) day 1.69193 0.05463 30.973 < 2e-16 *** mouseid01-01 7.13040 0.24045 29.655 < 2e-16 *** mouseid01-02 1.78189 0.28120 6.337 4.11e-07 *** mouseid01-03 6.49483 0.24045 27.011 < 2e-16 *** mouseid02-01 6.69657 0.24045 27.850 < 2e-16 *** mouseid02-02 5.84940 0.24045 24.327 < 2e-16 *** mouseid02-03 0.70050 0.33894 2.067 0.0469 * mouseid03-01 6.52969 0.24045 27.156 < 2e-16 *** mouseid03-02 6.09318 0.24045 25.341 < 2e-16 *** mouseid03-03 5.90790 0.24045 24.570 < 2e-16 *** mouseid04-01 6.58500 0.24045 27.386 < 2e-16 *** mouseid04-02 5.86158 0.24045 24.378 < 2e-16 *** mouseid04-03 6.32530 0.24045 26.306 < 2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.3958 on 32 degrees of freedom (30 observations deleted due to missingness) Multiple R-squared: 0.999, Adjusted R-squared: 0.9986 F-statistic: 2387 on 13 and 32 DF, p-value: < 2.2e-16 > > flow %>% + filter(day<=4) %>% + lm(log(Pd)~box:day+mouseid-1,data=.) -> fit2 > summary(fit2) Call: lm(formula = log(Pd) ~ box:day + mouseid - 1, data = .) Residuals: Min 1Q Median 3Q Max -0.54582 -0.21023 -0.01221 0.12066 0.49245 Coefficients: Estimate Std. Error t value Pr(>|t|) mouseid01-01 6.67620 0.29204 22.860 < 2e-16 *** mouseid01-02 1.23685 0.34621 3.572 0.00126 ** mouseid01-03 6.04063 0.29204 20.684 < 2e-16 *** mouseid02-01 6.40958 0.30581 20.959 < 2e-16 *** mouseid02-02 5.56240 0.30581 18.189 < 2e-16 *** mouseid02-03 0.29871 0.42943 0.696 0.49222 mouseid03-01 6.38173 0.27159 23.497 < 2e-16 *** mouseid03-02 5.94521 0.27159 21.890 < 2e-16 *** mouseid03-03 5.75993 0.27159 21.208 < 2e-16 *** mouseid04-01 7.29722 0.27159 26.868 < 2e-16 *** mouseid04-02 6.57380 0.27159 24.205 < 2e-16 *** mouseid04-03 7.03752 0.27159 25.912 < 2e-16 *** box01:day 1.87361 0.09602 19.512 < 2e-16 *** box02:day 1.80673 0.10265 17.600 < 2e-16 *** box03:day 1.75112 0.08589 20.389 < 2e-16 *** box04:day 1.40705 0.08589 16.383 3.38e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.3326 on 29 degrees of freedom (30 observations deleted due to missingness) Multiple R-squared: 0.9993, Adjusted R-squared: 0.999 F-statistic: 2747 on 16 and 29 DF, p-value: < 2.2e-16 > > anova(fit1,fit2) Analysis of Variance Table Model 1: log(Pd) ~ day + mouseid - 1 Model 2: log(Pd) ~ box:day + mouseid - 1 Res.Df RSS Df Sum of Sq F Pr(>F) 1 32 5.0131 2 29 3.2087 3 1.8044 5.4361 0.004321 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > > > ## ----beta-dose2---------------------------------------------------------- > coef(fit2) %>% + bind_rows() %>% + gather(var,val) %>% + mutate( + var=stri_replace_all_regex(var,"mouseid(\\d{2})-(\\d{2})","dose[$1-$2]"), + var=stri_replace_all_regex(var,"box(\\d{2}):day","Beta[$1]"), + val=exp(val) + ) -> theta1 > > expand.grid( + box=sprintf("%02d",1:5), + mouse=sprintf("%02d",1:3) + ) %>% + mutate( + mouseid=paste0(box,"-",mouse), + betavar=sprintf("Beta[%s]",box), + dosevar=sprintf("dose[%s]",mouseid) + ) %>% + left_join(theta1,by=c("betavar"="var")) %>% + rename(Beta=val) %>% + left_join(theta1,by=c("dosevar"="var")) %>% + rename(dose=val) %>% + select(-betavar,-dosevar) %>% + arrange(box,mouse) %>% + mutate( + Beta=coalesce(Beta,0), + dose=coalesce(dose,0) + ) -> theta > > theta box mouse mouseid Beta dose 1 01 01 01-01 6.511778 793.298642 2 01 02 01-02 6.511778 3.444740 3 01 03 01-03 6.511778 420.158396 4 02 01 02-01 6.090502 607.635623 5 02 02 02-02 6.090502 260.448444 6 02 03 02-03 6.090502 1.348118 7 03 01 03-01 5.761041 590.948395 8 03 02 03-02 5.761041 381.920847 9 03 03 03-03 5.761041 317.326831 10 04 01 04-01 4.083874 1476.190521 11 04 02 04-02 4.083874 716.085785 12 04 03 04-03 4.083874 1138.554819 13 05 01 05-01 0.000000 0.000000 14 05 02 05-02 0.000000 0.000000 15 05 03 05-03 0.000000 0.000000 > > > ## ----build_pomp,cache=FALSE,purl=TRUE------------------------------------ > library(doParallel) Loading required package: parallel > registerDoParallel() > > theta %>% + mutate( + sigmaPd = 2, + sigmaRBC = 0.1, + sigmaRetic = 0.3, + sigmaW = 1, + sigmaN = 0.5, + sigmaR = 0.5, + E_0 = 8.0e6, + R_0 = 3e5, + W_0 = 8.8e4, + N_0 = 8e5 + ) -> theta > > foreach (m = iter(theta,"row"),.inorder=TRUE,.combine=c) %dopar% { + + flow %>% + filter(mouseid==m$mouseid) %>% + select(day,Pd,RBC,Retic) %>% + mutate(Retic=if_else(day %in% c(0,14),NA_real_,Retic)) %>% + pomp( + params=select(m,-mouseid,-box,-mouse) %>% unlist(), + times="day", + t0=0, + rmeasure=Csnippet(" + Retic = rlnorm(log(1+R),sigmaRetic)-1; + RBC = rlnorm(log(1+E+R),sigmaRBC)-1; + Pd = rlnorm(log(1+K),sigmaPd)-1;"), + dmeasure=Csnippet(" + double l1, l2, l3; + l1 = (R_FINITE(Retic)) ? dlnorm(1+Retic,log(1+R),sigmaRetic,1) : 0; + l2 = (R_FINITE(RBC)) ? dlnorm(1+RBC,log(1+E+R),sigmaRBC,1) : 0; + l3 = (R_FINITE(Pd) && Pd>0) ? dlnorm(1+Pd,log(1+K),sigmaPd,1) : 0; + lik = (give_log) ? l1+l2+l3 : exp(l1+l2+l3);"), + rprocess=discrete_time( + step.fun=Csnippet(" + double Mold = M; + M = Beta*K*exp(-(W+N)/(R+E)); + E = (R+E)*exp(-(Mold+N)/(R+E)); + N = rlnorm(log(N),sigmaN); + W = rlnorm(log(W),sigmaW); + R = rlnorm(log(R),sigmaR); + K = (R+E>0) ? (R+E)*(1-exp(-M/(R+E))): 0; + "), + delta.t=1 + ), + partrans=parameter_trans( + log=c("sigmaW","sigmaR","sigmaN", + "sigmaPd","sigmaRBC","sigmaRetic", + "N_0","W_0","E_0","R_0") + ), + rinit=Csnippet(" + E = E_0; + R = R_0; + N = N_0; + W = W_0; + M = 0; + K = dose;"), + statenames=c("E","R","W","N","M","K"), + paramnames=c( + "Beta","dose", + "sigmaPd","sigmaRBC","sigmaRetic", + "sigmaW","sigmaN","sigmaR", + "E_0","R_0","W_0","N_0" + ) + ) + } %>% + set_names(theta$mouseid) -> pos > > > ## ----mpi_setup,purl=TRUE,echo=FALSE-------------------------------------- > if (file.exists("CLUSTER")) { + scan("CLUSTER",what=integer(0)) -> ncpu + library(doMPI) + cl <- startMPIcluster(ncpu,verbose=TRUE,logdir="/tmp") + registerDoMPI(cl) + } Read 1 item Loading required package: Rmpi Master processor name: theorygroup01; nodename: theorygroup01 Size of MPI universe: 1 Spawning 250 workers using the command: /export/local/apps/R/R-3.6.1/lib64/R/bin/Rscript /userdata/kingaa/R/x86_64-pc-linux-gnu-library/3.6/doMPI/RMPIworker.R WORKDIR=/home/kingaa/code LOGDIR=/tmp MAXCORES=1 COMM=3 INTERCOMM=4 MTAG=10 WTAG=11 INCLUDEMASTER=TRUE BCAST=TRUE VERBOSE=TRUE App launch reported: 12 (out of 12) daemons - 229 (out of 250) procs 250 slaves are spawned successfully. 0 failed. > > > ## ----shared_params,purl=TRUE--------------------------------------------- > shared_params <- c( + "E_0","N_0","R_0", + "sigmaN","sigmaPd","sigmaR", + "sigmaRBC","sigmaRetic","sigmaW", + "W_0" + ) > > > ## ----mf1,cache=FALSE,purl=TRUE------------------------------------------- > bake(file="m5mf1.rds",{ + + theta %>% + select(-box,-mouse,-mouseid,-Beta,-dose) %>% + gather(var,val) %>% + group_by(var) %>% + summarize(min=min(val),max=max(val)) %>% + ungroup() %>% + gather(bound,val,-var) %>% + spread(var,val) %>% + column_to_rownames("bound") %>% + as.matrix() -> box + + profileDesign( + mouseid=theta$mouseid, + lower=0.5*box["min",shared_params], + upper=1.5*box["max",shared_params], + nprof=50 + ) -> starts + + registerDoRNG(2022051968) + + foreach (start=iter(starts,"row"), + .inorder=FALSE,.combine=bind_rows, + .packages=c("pomp","magrittr","plyr","dplyr") + ) %dopar% + { + tryCatch({ + pos[[start$mouseid]] -> m + start %>% select(-mouseid) %>% unlist() -> theta + coef(m,names(theta)) <- unname(theta) + start <- cbind(start["mouseid"],as.data.frame(as.list(coef(m)))) + m %>% + mif2(Nmif=50,Np=10000,tol=exp(-60), + cooling.fraction.50=0.1,cooling.type="geometric", + rw.sd=rw.sd( + sigmaPd=0.02,sigmaRBC=0.02,sigmaRetic=0.02, + sigmaW=0.02,sigmaN=0.02,sigmaR=0.02, + E_0=ivp(0.05),R_0=ivp(0.05),W_0=ivp(0.05),N_0=ivp(0.05) + ) + ) %>% + mif2(Nmif=100) -> m + raply(5,m %>% pfilter(Np=50000,tol=exp(-60)) %>% logLik()) %>% + logmeanexp(se=TRUE) -> pfs + m %>% coef() %>% rbind() %>% as.data.frame() -> fin + fin$loglik <- pfs[[1]] + fin$loglik.se <- pfs[[2]] + fin$mouseid <- start$mouseid + fin$msg <- "mf1" + start$msg <- "start" + bind_rows(start,fin) + }, + error=function (e) { + fin <- start + fin$msg <- conditionMessage(e) + full_join(start,fin,by=names(start)) + }) + } -> res + attr(res,"ncpu") <- getDoParWorkers() + + res + }) -> mf1 > > > > > > > ## ----exclude------------------------------------------------------------- > control_mice <- grep("^05-",names(pos),perl=TRUE) > > > ## ----mf2,cache=FALSE,purl=TRUE------------------------------------------- > bake(file="m5mf2.rds",{ + + mf1 %>% + separate(mouseid,into=c("box","mouse"),remove=FALSE) %>% + ## remove controls and mouse that died early + filter(box != "05" & mouseid != "01-01") %>% + group_by(mouseid) %>% + filter(!is.na(loglik) | loglik>max(loglik,na.rm=TRUE)-10) %>% + ungroup() %>% + select(starts_with("sigma"),ends_with("_0")) %>% + gather(var,val) %>% + group_by(var) %>% + summarize(min=min(val),max=max(val)) %>% + ungroup() %>% + gather(bound,val,-var) %>% + spread(var,val) %>% + column_to_rownames("bound") %>% + as.matrix() -> box + + sobolDesign( + lower=box["min",shared_params], + upper=box["max",shared_params], + nseq=250 + ) -> starts + + registerDoRNG(1510289681) + + foreach (start=iter(starts,"row"), + .inorder=FALSE,.combine=bind_rows, + .packages=c("panelPomp","magrittr","plyr","dplyr"), + .options.mpi=list(chunkSize=1)) %dopar% + { + tryCatch({ + panelPomp(pos[-control_mice],shared=unlist(start)) -> m + m %>% coef() %>% rbind() %>% as.data.frame() -> start + m %>% + mif2(Nmif=300,Np=20000,tol=exp(-180), + cooling.fraction.50=0.75,cooling.type="geometric", + rw.sd=rw.sd( + sigmaPd=0.01,sigmaRBC=0.01,sigmaRetic=0.01, + sigmaW=0.01,sigmaN=0.01,sigmaR=0.01, + E_0=ivp(0.02),R_0=ivp(0.02),W_0=ivp(0.02),N_0=ivp(0.02) + ) + ) -> m + m %>% coef() %>% rbind() %>% as.data.frame() -> fin + fin$msg <- "mf2" + start$msg <- "mf1" + bind_rows(start,fin) + }, + error=function (e) { + fin <- start + fin$msg <- conditionMessage(e) + full_join(start,fin,by=names(start)) + }) + } -> res + attr(res,"ncpu") <- getDoParWorkers() + + res + }) -> mf2 > > > > > ## ----mf3,cache=FALSE,purl=TRUE------------------------------------------- > bake(file="m5mf3.rds",{ + + mf2 %>% + filter(msg=="mf2") %>% + arrange(sigmaW,sigmaN,sigmaR) %>% + select(-msg) -> starts + + registerDoRNG(599108918) + + foreach (start=iter(starts,"row"), + .inorder=FALSE,.combine=bind_rows, + .packages=c("panelPomp","magrittr","plyr","dplyr"), + .options.mpi=list(chunkSize=1)) %dopar% + { + tryCatch({ + panelPomp(pos[-control_mice],shared=unlist(start)) -> m + m %>% coef() %>% rbind() %>% as.data.frame() -> start + m %>% + mif2(Nmif=100,Np=20000,tol=exp(-180), + cooling.fraction.50=0.1,cooling.type="geometric", + rw.sd=rw.sd( + sigmaPd=0.01,sigmaRBC=0.01,sigmaRetic=0.01, + sigmaW=0.01,sigmaN=0.01,sigmaR=0.01, + E_0=ivp(0.02),R_0=ivp(0.02),W_0=ivp(0.02),N_0=ivp(0.02) + ) + ) -> m + m %>% coef() %>% rbind() %>% as.data.frame() -> fin + fin$msg <- "mf3" + start$msg <- "mf2" + bind_rows(start,fin) + }, + error=function (e) { + fin <- start + fin$msg <- conditionMessage(e) + full_join(start,fin,by=names(start)) + }) + } -> res + attr(res,"ncpu") <- getDoParWorkers() + + res + }) -> mf3 > > > > > ## ----pf3,cache=FALSE,purl=TRUE------------------------------------------- > bake(file="m5pf3.rds",{ + + bind_rows(mf2,mf3) %>% + distinct() %>% + arrange(sigmaW,sigmaN,sigmaR) -> starts + + registerDoRNG(865465456) + + foreach (start=iter(starts,"row"), + .inorder=FALSE,.combine=bind_rows, + .packages=c("panelPomp","magrittr","plyr","dplyr"), + .options.mpi=list(chunkSize=1)) %dopar% + { + tryCatch({ + start %>% select(-msg) %>% unlist() -> p + panelPomp(pos[-control_mice],shared=p) -> m + replicate( + n=10, + m %>% pfilter(Np=150000,tol=exp(-180)) + ) %>% + sapply(logLik) %>% + logmeanexp(se=TRUE) -> lls + start$loglik <- lls[1] + start$loglik.se <- lls[2] + start + }, + error=function (e) { + start$msg <- conditionMessage(e) + start + }) + } -> res + attr(res,"ncpu") <- getDoParWorkers() + + res + }) -> pf3 > > > > > ## ----pf3-2--------------------------------------------------------------- > pf3 %>% + group_by(msg) %>% + filter(loglik+2*loglik.se>max(loglik-2*loglik.se)-100) %>% + summarize(n=length(loglik),max=max(loglik), + mean=mean(loglik),sd=sd(loglik),meanSE=mean(loglik.se)) # A tibble: 3 x 6 msg n max mean sd meanSE 1 mf1 42 -11312. -11388. 23.0 3.18 2 mf2 246 -11227. -11271. 23.0 4.38 3 mf3 240 -11225. -11268. 22.9 4.28 > > > ## ----mf4,purl=TRUE,cache=FALSE------------------------------------------- > bake(file="m5mf4.rds",{ + + shared_params <- c( + "E_0","N_0","R_0", + "sigmaN","sigmaPd","sigmaR", + "sigmaRBC","sigmaRetic", + "W_0" + ) + + pf3 %>% + as_tibble() %>% + filter(loglik>max(loglik)-100) %>% + select(starts_with("sigma"),ends_with("_0")) %>% + gather(var,val) %>% + group_by(var) %>% + summarize(min=min(val),max=max(val)) %>% + ungroup() %>% + gather(bound,val,-var) %>% + spread(var,val) %>% + column_to_rownames("bound") %>% + as.matrix() -> box + + profileDesign( + sigmaW=seq(from=0.2,to=box["max","sigmaW"],length=100), + lower=box["min",shared_params], + upper=box["max",shared_params], + nprof=10 + ) -> starts + + registerDoRNG(933950321) + + foreach ( + start=iter(starts,"row"), + .combine=bind_rows, + .packages=c("plyr","dplyr","magrittr","panelPomp") + ) %dopar% { + + start %>% unlist() -> p + panelPomp(pos[-control_mice],shared=p) -> m + m %>% + mif2(Nmif=100,Np=10000,tol=exp(-180), + rw.sd=rw.sd( + sigmaPd=0.02,sigmaRBC=0.02,sigmaRetic=0.02, + sigmaW=0,sigmaN=0.02,sigmaR=0.02, + E_0=ivp(0.02),R_0=ivp(0.02),W_0=ivp(0.02),N_0=ivp(0.02) + ), + cooling.type="geometric",cooling.fraction.50=0.1) %>% + mif2(Np=20000) %>% + mif2(Np=40000) -> m + m %>% coef() %>% rbind() %>% as.data.frame() %>% + cbind(data.frame(msg="mf4")) + } -> res + attr(res,"ncpu") <- getDoParWorkers() + res + + }) -> mf4 > > > ## ----pf4,purl=TRUE,cache=FALSE------------------------------------------- > bake(file="m5pf4.rds",{ + + mf4 %>% + arrange(sigmaW,sigmaN,sigmaR) %>% + mutate(id=row_number(msg)) -> starts + + registerDoRNG(537031349) + + foreach ( + start=iter(starts,"row"), + .combine=bind_rows, + .packages=c("plyr","dplyr","tibble","magrittr","panelPomp") + ) %dopar% { + + start %>% select(-msg) %>% unlist() -> p + panelPomp(pos[-control_mice],shared=p) -> m + replicate(n=10, + m %>% pfilter(Np=150000,tol=exp(-180)) %>% logLik() + ) -> lls + start %>% + remove_rownames() %>% + cbind(data.frame(loglik=lls)) + + } -> res + attr(res,"ncpu") <- getDoParWorkers() + res + + }) -> pf4 > > > > > ## ----mle----------------------------------------------------------------- > pf4 %>% + select(loglik,starts_with("sigma"),ends_with("_0")) %>% + filter(loglik==max(loglik)) -> mle > > mle %>% + gather(var,val,-loglik) %>% + mutate(val=signif(val,3)) %>% + spread(var,val) %>% + select(loglik,starts_with("sigma"),ends_with("_0")) loglik sigmaN sigmaPd sigmaR sigmaRBC sigmaRetic sigmaW E_0 N_0 1 -11217.79 0.625 0.597 0.393 0.0794 0.151 1.04 7720000 411000 R_0 W_0 1 618000 176000 > > > ## ----sm1,cache=FALSE,purl=TRUE------------------------------------------- > bake(file="m5sm1.rds",{ + + registerDoRNG(1510516029) + + foreach (i=1:2000, + .inorder=FALSE,.combine=bind_rows, + .packages=c("panelPomp","magrittr","reshape2","plyr","tibble","tidyr","dplyr") + ) %dopar% + { + + mle %>% select(-loglik) %>% unlist() -> p + panelPomp(pos,shared=p) -> m + + m %>% + pfilter(Np=1e5,filter.traj=TRUE,tol=exp(-180)) %>% + as("list") %>% + llply(filter.traj) %>% + melt() %>% + rename(mouseid=L1) %>% + mutate(rep=i) + + } -> res + attr(res,"ncpu") <- getDoParWorkers() + res + + }) -> sm1 > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > ## ----session-info,cache=FALSE-------------------------------------------- > bake(file="session_info.rds",{ + sessionInfo() + }) %>% + print() R version 3.6.1 (2019-07-05) Platform: x86_64-pc-linux-gnu (64-bit) Running under: Ubuntu 18.04.3 LTS Matrix products: default BLAS: /usr/lib/x86_64-linux-gnu/openblas/libblas.so.3 LAPACK: /usr/lib/x86_64-linux-gnu/libopenblasp-r0.2.20.so locale: [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8 [4] LC_COLLATE=C LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 [7] LC_PAPER=en_US.UTF-8 LC_NAME=C LC_ADDRESS=C [10] LC_TELEPHONE=C LC_MEASUREMENT=C LC_IDENTIFICATION=C attached base packages: [1] parallel stats graphics grDevices utils datasets methods [8] base other attached packages: [1] doMPI_0.2.2 Rmpi_0.6-9 doParallel_1.0.15 aakmisc_0.27-1 [5] doRNG_1.7.1 rngtools_1.4 pkgmaker_0.27 registry_0.5-1 [9] iterators_1.0.12 foreach_1.4.7 panelPomp_0.9.1.3 pomp_2.3.0.5 [13] stringi_1.4.3 magrittr_1.5 forcats_0.4.0 stringr_1.4.0 [17] dplyr_0.8.3 purrr_0.3.2 readr_1.3.1 tidyr_1.0.0 [21] tibble_2.1.3 ggplot2_3.2.1 tidyverse_1.2.1 plyr_1.8.4 loaded via a namespace (and not attached): [1] Rcpp_1.0.2 lubridate_1.7.4 mvtnorm_1.0-11 lattice_0.20-38 [5] utf8_1.1.4 assertthat_0.2.1 zeallot_0.1.0 digest_0.6.21 [9] R6_2.4.0 cellranger_1.1.0 backports_1.1.5 coda_0.19-3 [13] httr_1.4.1 pillar_1.4.2 rlang_0.4.0 lazyeval_0.2.2 [17] RPostgreSQL_0.6-2 curl_4.2 readxl_1.3.1 rstudioapi_0.10 [21] munsell_0.5.0 broom_0.5.2 compiler_3.6.1 modelr_0.1.5 [25] pkgconfig_2.0.3 tidyselect_0.2.5 codetools_0.2-16 fansi_0.4.0 [29] crayon_1.3.4 withr_2.1.2 grid_3.6.1 nlme_3.1-141 [33] jsonlite_1.6 xtable_1.8-4 gtable_0.3.0 lifecycle_0.1.0 [37] DBI_1.0.0 scales_1.0.0 bibtex_0.4.2 cli_1.1.0 [41] reshape2_1.4.3 xml2_1.2.2 ellipsis_0.3.0 generics_0.0.2 [45] vctrs_0.2.0 deSolve_1.24 tools_3.6.1 glue_1.3.1 [49] hms_0.5.1 colorspace_1.4-1 rvest_0.3.4 haven_2.1.1 > > > proc.time() user system elapsed 2586.18 34039.71 36617.48