## ----style, echo = FALSE, results = 'asis'------------------------------------ # BiocStyle::markdown() ## ----include = FALSE---------------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>") ## ----setup, eval=TRUE, echo=FALSE, include=FALSE----------------------------- library(bmstdr) library(ggplot2) library(tidyr) library(RColorBrewer) knitr::opts_chunk$set(eval = F) longrun <- FALSE tabs <- lapply(list.files(system.file('txttables', package = 'bmstdr'), full.names = TRUE), dget) # fns <- list.files(system.file('txttables', package = 'bmstdr')) # how the table data and file names match table1 <- tabs[[1]] table2 <- tabs[[4]] table3 <- tabs[[5]] table4 <- tabs[[6]] table5 <- tabs[[7]] table6 <- tabs[[8]] table7 <- tabs[[9]] table8 <- tabs[[10]] table9 <- tabs[[11]] table10 <- tabs[[2]] table11 <- tabs[[3]] tablepath <- system.file('last3tables', package = 'bmstdr') print(tablepath) table4.1 <- dget(file=paste0(tablepath, "/table4.1.txt")) table4.2 <- dget(file=paste0(tablepath, "/table4.2.txt")) table4.3 <- dget(file=paste0(tablepath, "/table4.3.txt")) ## ----------------------------------------------------------------------------- # M1 <- Bspatial(formula=yo3~xmaxtemp+xwdsp+xrh, data=nyspatial, mchoice=T) ## ----------------------------------------------------------------------------- # M2 <- Bspatial(model="spat", formula=yo3~xmaxtemp+xwdsp+xrh, data=nyspatial, # coordtype="utm", coords=4:5, phi=0.4, mchoice=T) ## ----echo=TRUE, eval=FALSE---------------------------------------------------- # asave <- phichoice_sp() # print(asave) ## ----echo=TRUE, eval=FALSE---------------------------------------------------- # M3 <- Bspatial(package="spBayes", formula=yo3~xmaxtemp+xwdsp+xrh, data=nyspatial, # coordtype="utm", coords=4:5, prior.phi=c(0.005, 2), mchoice=T) # M4 <- Bspatial(package="stan", formula=yo3~xmaxtemp+xwdsp+xrh, data=nyspatial, # coordtype="utm", coords=4:5,phi=0.4, mchoice=T) # M5 <- Bspatial(package="inla",formula=yo3~xmaxtemp+xwdsp+xrh, data=nyspatial, # coordtype="utm", coords=4:5, mchoice=T) ## ----echo=FALSE, eval=FALSE--------------------------------------------------- # a3 <- Bmchoice(case="MC.sigma2.unknown", y=ydata) # # Now organize the all the results for forming Table 1. # a5 <- rep(NA, 11) # a5[c(1, 3, 5, 7, 9:11)] <- unlist(M5$mchoice) # table1 <- cbind.data.frame(unlist(a3), M1$mchoice, M2$mchoice, M3$mchoice, M4$mchoice, a5) # colnames(table1) <- paste("M", 0:5, sep="") ## ----echo=TRUE, eval=FALSE---------------------------------------------------- # Bmchoice(case="MC.sigma2.unknown", y=ydata). ## ----echo=TRUE, eval=FALSE---------------------------------------------------- # s <- c(8,11,12,14,18,21,24,28) # f1 <- yo3~xmaxtemp+xwdsp+xrh # M1.v <- Bspatial(package="none", model="lm", formula=f1, # data=nyspatial, validrows=s) # M2.v <- Bspatial(package="none", model="spat", formula=f1, # data=nyspatial, coordtype="utm", coords=4:5,phi=0.4, validrows=s) # M3.v <- Bspatial(package="spBayes", prior.phi=c(0.005, 2), formula=f1, # data=nyspatial, coordtype="utm", coords=4:5, validrows=s) # M4.v <- Bspatial(package="stan",formula=f1, # data=nyspatial, coordtype="utm", coords=4:5,phi=0.4 , validrows=s) # M5.v <- Bspatial(package="inla", formula=f1, data=nyspatial, # coordtype="utm", coords=4:5, validrows=s) ## ----echo=TRUE, eval=TRUE----------------------------------------------------- set.seed(44) x <- runif(n=28) u <- order(x) s1 <- u[1:7] s2 <- u[8:14] s3 <- u[15:21] s4 <- u[22:28] ## ----valplot1, echo=T, eval=FALSE, message=FALSE, results='hide', fig.cap="Prediction against observation plot with the prediction intervals included. The `in/out' symbol in the plot indicates whether or not a prediction interval incudes the 45 degree line."---- # M2.v3 <- Bspatial(model="spat", formula=yo3~xmaxtemp+xwdsp+xrh, data=nyspatial, # coordtype="utm", coords=4:5, validrows= s3, phi=0.4, verbose = FALSE) ## ----echo=T, eval=FALSE, message=FALSE, results='hide'------------------------ # names(M2.v3) # psums <- get_validation_summaries(M2.v3$valpreds) # names(psums) # a <- obs_v_pred_plot(yobs=M2.v3$yobs_preds$yo3, predsums=psums, segments=F, summarystat = "mean" ) ## ----echo=TRUE, eval=FALSE, message=FALSE, results='hide'--------------------- # f2 <- y8hrmax~xmaxtemp+xwdsp+xrh # M1 <- Bsptime(model="lm", formula=f2, data=nysptime, scale.transform = "SQRT", N=5000) # M2 <- Bsptime(model="separable", formula=f2, data=nysptime, scale.transform = "SQRT", # coordtype="utm", coords=4:5, N=5000) ## ----echo=TRUE, eval=FALSE, message=FALSE, results='hide'--------------------- # a <- residuals(M1, numbers=list(sn=28, tn=62)) ## ----ptime, echo=T, eval=T, fig.cap="Weekly Covid-19 death rate per 100,000"---- engdeaths$covidrate <- 100000*engdeaths$covid/engdeaths$popn ptime <- ggplot(data=engdeaths, aes(x=factor(Weeknumber), y=covidrate)) + geom_boxplot() + labs(x = "Week", y = "Death rate per 100,000") + stat_summary(fun=median, geom="line", aes(group=1, col="red")) + theme(legend.position = "none") ptime ## ----echo=T, eval=T, message=FALSE, results='hide'---------------------------- Ncar <- 50000 burn.in.car <- 10000 thin <- 10 ## ----echo=TRUE, eval=FALSE, message=FALSE, results='hide'--------------------- # f1 <- noofhighweeks ~ jsa + log10(houseprice) + log(popdensity) + sqrt(no2) ## ----echo=T, eval=FALSE, message=FALSE, results='hide'------------------------ # M1 <- Bcartime(formula=f1, data=engtotals, family="binomial", # trials=engtotals$nweek, N=Ncar, burn.in=burn.in.car, thin=thin) ## ----echo=T, eval=FALSE, message=FALSE, results='hide'------------------------ # M1.leroux <- Bcartime(formula=f1, data=engtotals, scol="spaceid", # model="leroux", W=Weng, family="binomial", trials=engtotals$nweek, # N=Ncar, burn.in=burn.in.car, thin=thin) ## ----echo=T, eval=FALSE, message=FALSE, results='hide'------------------------ # M1.bym <- Bcartime(formula=f1, data=engtotals, # scol="spaceid", model="bym", W=Weng, family="binomial", # trials=engtotals$nweek, N=Ncar, burn.in=burn.in.car, thin=thin) ## ----echo=T, eval=FALSE, message=FALSE, results='hide'------------------------ # M1.inla.bym <- Bcartime(formula=f1, data=engtotals, scol ="spaceid", # model=c("bym"), W=Weng, family="binomial", trials=engtotals$nweek, # package="inla", N=Ncar, burn.in=burn.in.car, thin=thin) ## ----echo=T, eval=FALSE, message=FALSE, results='hide'------------------------ # a <- rbind(M1$mchoice, M1.leroux$mchoice, M1.bym$mchoice) # a <- a[, -(5:6)] # a <- a[, c(2, 1, 4, 3)] # b <- M1.inla.bym$mchoice[1:4] # a <- rbind(a, b) # rownames(a) <- c("Independent", "Leroux", "BYM", "INLA-BYM") # colnames(a) <- c("pDIC", "DIC", "pWAIC", "WAIC") # table4.1 <- a ## ----echo=T, eval=FALSE------------------------------------------------------- # f2 <- covid ~ offset(logEdeaths) + jsa + log10(houseprice) + log(popdensity) + sqrt(no2) ## ----echo=T, eval=FALSE, message=FALSE, results='hide'------------------------ # M2 <- Bcartime(formula=f2, data=engtotals, family="poisson", # N=Ncar, burn.in=burn.in.car, thin=thin) # # M2.leroux <- Bcartime(formula=f2, data=engtotals, # scol="spaceid", model="leroux", family="poisson", W=Weng, # N=Ncar, burn.in=burn.in.car, thin=thin) # # M2.bym <- Bcartime(formula=f2, data=engtotals, # scol="spaceid", model="bym", family="poisson", W=Weng, # N=Ncar, burn.in=burn.in.car, thin=thin) # # M2.inla.bym <- Bcartime(formula=f2, data=engtotals, scol ="spaceid", # model=c("bym"), family="poisson", # W=Weng, offsetcol="logEdeaths", link="log", # package="inla", N=Ncar, burn.in = burn.in.car, thin=thin) # ## ----echo=FALSE, eval=FALSE, message=FALSE, results='hide'-------------------- # a <- rbind(M2$mchoice, M2.leroux$mchoice, M2.bym$mchoice) # a <- a[, -(5:6)] # a <- a[, c(2, 1, 4, 3)] # b <- M2.inla.bym$mchoice[1:4] # a <- rbind(a, b) # rownames(a) <- c("Independent", "Leroux", "BYM", "INLA-BYM") # colnames(a) <- c("pDIC", "DIC", "pWAIC", "WAIC") # table4.2 <- a # dput(table4.2, file=paste0(tablepath, "/table4.2.txt")) ## ----echo=T, eval=FALSE------------------------------------------------------- # f3 <- sqrt(no2) ~ jsa + log10(houseprice) + log(popdensity) ## ----echo=T, eval=FALSE, message=FALSE, results='hide'------------------------ # M3 <- Bcartime(formula=f3, data=engtotals, family="gaussian", # N=Ncar, burn.in=burn.in.car, thin=thin) # # M3.leroux <- Bcartime(formula=f3, data=engtotals, # scol="spaceid", model="leroux", family="gaussian", W=Weng, # N=Ncar, burn.in=burn.in.car, thin=thin) # # # M3.inla.bym <- Bcartime(formula=f3, data=engtotals, scol ="spaceid", # model=c("bym"), family="gaussian", # W=Weng, package="inla", N=Ncar, burn.in =burn.in.car, thin=thin) # ## ----echo=T, eval=FALSE, message=FALSE, results='hide'------------------------ # a <- rbind(M3$mchoice, M3.leroux$mchoice) # a <- a[, -(5:6)] # a <- a[, c(2, 1, 4, 3)] # b <- M3.inla.bym$mchoice[1:4] # a <- rbind(a, b) # rownames(a) <- c("Independent", "Leroux", "INLA-BYM") # colnames(a) <- c("pDIC", "DIC", "pWAIC", "WAIC") # table4.3 <- a # dput(table4.3, file=paste0(tablepath, "/table4.3.txt")) ## ----echo=TRUE, eval=FALSE---------------------------------------------------- # nweek <- rep(1, nrow(engdeaths)) ## ----echo=TRUE, eval=FALSE---------------------------------------------------- # f1 <- highdeathsmr ~ jsa + log10(houseprice) + log(popdensity) ## ----echo=T, eval=FALSE------------------------------------------------------- # M1st <- Bcartime(formula=f1, data=engdeaths, scol=scol, tcol=tcol, trials=nweek, # W=Weng, model="linear", family="binomial", package="CARBayesST", N=Ncar, # burn.in=burn.in.car, thin=thin) ## ----echo=TRUE, eval=FALSE---------------------------------------------------- # f2 <- covid ~ offset(logEdeaths) + jsa + log10(houseprice) + log(popdensity) + n0 + n1 + n2 + n3 ## ----echo=FALSE, eval=FALSE--------------------------------------------------- # M2st <- Bcartime(formula=f2, data=engdeaths, scol=scol, tcol=tcol, W=Weng, # model="ar", family="poisson", package="CARBayesST", N=Ncar, burn.in=burn.in.car, thin=thin) ## ----echo=TRUE, eval=FALSE---------------------------------------------------- # vs <- sample(nrow(engdeaths), 0.1*nrow(engdeaths)) ## ----obsvpredplot-m2st, echo=TRUE, eval=FALSE, message=FALSE, results='hide', fig.show='hide'---- # f20 <- covid ~ offset(logEdeaths) + jsa + log10(houseprice) + log(popdensity) + n0 # model <- c("bym", "ar1") # f2inla <- covid ~ jsa + log10(houseprice) + log(popdensity) + n0 # set.seed(5) # vs <- sample(nrow(engdeaths), 0.1*nrow(engdeaths)) # M2st_ar2.0 <- Bcartime(formula=f20, data=engdeaths, scol="spaceid", tcol= "Weeknumber", # W=Weng, model="ar", AR=2, family="poisson", package="CARBayesST", # N=Ncar, burn.in=burn.in.car, thin=thin, # validrows=vs, verbose=F) # M2stinla.0 <- Bcartime(data=engdeaths, formula=f2inla, W=Weng, scol ="spaceid", tcol="Weeknumber", # offsetcol="logEdeaths", model=model, link="log", family="poisson", package="inla", validrow=vs, N=N, burn.in=0) # yobspred <- M2st_ar2.0$yobs_preds # names(yobspred) # yobs <- yobspred$covid # predsums <- get_validation_summaries(t(M2st_ar2.0$valpreds)) # dim(predsums) # b <- obs_v_pred_plot(yobs, predsums, segments=T) # names(M2stinla.0) # inlapredsums <- get_validation_summaries(t(M2stinla.0$valpreds)) # dim(inlapredsums) # a <- obs_v_pred_plot(yobs, inlapredsums, segments=T) # inlavalid <- a$pwithseg # ar2valid <- b$pwithseg # library(ggpubr) # ggarrange(ar2valid, inlavalid, common.legend = TRUE, legend = "top", nrow = 2, ncol = 1) # figpath <- system.file("figures", package = "bmstdr") # ggsave(filename = paste0(figpath, "/figure11.png")) ## ----obsvpredplot, echo=FALSE, eval=TRUE, fig.cap="Predictions with 95% limits against observations for two models: AR (2) on the left panel and INLA on the right panel", fig.width=1.2---- figpath <- system.file("figures", package = "bmstdr") knitr::include_graphics(paste0(figpath, "/figure11.png")) ## ----echo=T, eval=FALSE------------------------------------------------------- # M3st <- Bcartime(formula=f3, data=engdeaths, scol=scol, tcol=tcol, # W=Weng, model="ar", family="gaussian", package="CARBayesST", # N=Ncar, burn.in=burn.in.car, thin=thin)