inst/applications/ex_DVC.R

source("setup.R")
load("DVC.Rda")

DVC <- data.frame(number = c(rowSums(x06), rowSums(x09)),
                  days = c(1:365, 1:365),
                  year = factor(c(rep("2006", 365), rep("2009", 365))))

mod <- mctm(bbs(number, df = 3, differences = 1) ~ 
                bbs(days, df = 3, differences = 1, cyclic = TRUE)
              + bbs(days, by = year, differences = 1, df = 3, cyclic = TRUE),
              data = DVC, family = Binomial(link = "probit"), ngrid = 0:20 * 10,
              control = boost_control(nu = 0.5, mstop = 250))

folds <- cv(model.weights(mod), strata = DVC$year)
(cv <- cvrisk(mod, folds = folds))
mod[mstop(cv)]

p <- predict(mod, newdata = DVC, annotate = TRUE, type = "response")
DVC$ID <- 1:nrow(DVC)
DVC$days <- c(seq(as.Date("2006-01-01"), length.out=365, by="1 day"),
              seq(as.Date("2009-01-01"), length.out=365, by="1 day"))
p <- merge(p, DVC, by = "ID")

save(p, file = "ex_DVC.Rda")

Try the ctmDevel package in your browser

Any scripts or data that you put into this service are public.

ctmDevel documentation built on May 2, 2019, 4:52 p.m.