################################################################################# ## 2. Getting started ################################################################################# library("informR") data("atus80ord", package = "informR") data("atus80int", package = "informR") atus80ord[1:5, ] atus80int[1:10, ] atus80ord[which(is.na(atus80ord[, "Activities"])), "Activities"] <- "MISSING" rawevents <- cbind(atus80ord$Activities, atus80ord$TUCASEID) evls <- gen.evl(rawevents, null.events = c("MISSING")) names(evls) evls$eventlist[[1]] evls$event.key evls$null.events alpha.ints <- gen.intercepts(evls, basecat = "Sleeping") length(alpha.ints) dim(alpha.ints[[1]][[1]]) alpha.ints[[1]][[1]][1, , ] alpha.fit <- rem(eventlist = evls$eventlist, statslist = alpha.ints, estimator = "BPM", prior.param = list(mu = 0, sigma = 100, nu = 4)) summary(alpha.fit) pois.mle <- log(prop.table(table(atus80ord$Activities))[-c(7, 10)]/prop.table(table(atus80ord$Activities))[10]) round(cbind(BPM = alpha.fit$coef[order(names(alpha.fit$coef))], pois.mle), 4) ################################################################################# ## 3. Sequence statistics ################################################################################# a1 <- paste(evls$event.key[-9, 1], evls$event.key[-9, 1], sep = "") beta.sforms <- gen.sformlist(evls, a1) beta.sforms[[1]][1, , ] evls$eventlist$"20030101033341" sapply(attr(evls$eventlist$"20030101033341", "char")[1:3], sf2nms, event.key = evls$event.key) beta.sforms$"20030101033341"[1:4, , c("aa", "bb")] beta.ints <- slbind(beta.sforms, alpha.ints, new.names = TRUE, event.key = evls$event.key) beta.fit <- rem(evls$eventlist, beta.ints, estimator = "BPM", prior.param = list(mu = 0, sigma = 100, nu = 4))$coef round(cbind(beta.fit[13:25]), 4) comComHaz <- exp(sum(beta.fit[c(5, 18)])) allOthHaz <- exp(beta.fit[(1:12)[-5]]) comComHaz/sum(comComHaz + allOthHaz) ################################################################################# ## 3.1. Constructing complex sequences ################################################################################# a2 <- c("a(c|h)a", "al+a") gamma.sforms <- gen.sformlist(evls, a2) gamma.ints <- slbind(gamma.sforms, beta.ints, new.names = TRUE, event.key = evls$event.key) gamma.fit <- rem(evls$eventlist, gamma.ints, estimator = "BPM", prior.param = list(mu = 0, sigma = 100, nu = 4)) summary(gamma.fit) sleep.sfs <- paste("a", letters[2:14], "a", sep = "") sleep.sforms <- glb.sformlist(evls, sforms = list(sleep.sfs), new.names = "InterSleep") sleep.ints <- slbind(sleep.sforms, gamma.ints) sleep.fit <- rem(evls$eventlist, sleep.ints, estimator = "BPM", prior.param = list(mu = 0, sigma = 100, nu = 4)) round(cbind(BPM = sleep.fit$coef, Z = sleep.fit$coef/sleep.fit$sd)[26:28, ], 4) ################################################################################# ## 4. Additional s-form parametrizations ################################################################################# tmp.sforms <- glb.sformlist(evls, list(c("a(b+|c)d", "a(b|c+)d")), new.names = "a(b+|c+)d") evls$eventlist[[1]] tmp.sforms[[1]] ################################################################################# ## 4.1. Conditioning out the prefix ################################################################################# abc.sform <- gen.sformlist(evls, c("abc"), cond = FALSE) abc.cond.sform <- gen.sformlist(evls, c("abc"), cond = TRUE) abc.sform[[1]] abc.cond.sform[[1]] ################################################################################# ## 4.2. Modifying statslists ################################################################################# sforms <- c("bf", "ef", "gf", "cf", "egf", "cef", "af") delta.sforms <- gen.sformlist(evls, sforms) delta.ints <- slbind(delta.sforms, alpha.ints) delta.fit <- rem(evls$eventlist, delta.ints) summary(delta.fit) delta.ints2 <- sldrop(delta.ints, c("egf", "Travel")) delta.fit2 <- rem(evls$eventlist, delta.ints2) names(delta.fit2$coef) c(delta.fit$BIC, delta.fit2$BIC) sl.ind <- 26:27 fem.evs <- unlist(glapply(atus80ord$SEX, atus80ord$TUCASEID, unique, regroup = FALSE)) fem.evs <- ifelse(fem.evs == 1, 1, 0) gamma.ints2 <- slbind.cond(fem.evs, gamma.ints, sl.ind = sl.ind, var.suffix = "FEM") sl.ind <- 26:29 gamma.fit2 <- rem(evls$eventlist, gamma.ints2, estimator = "BPM", prior.param = list(mu = 0, sigma = 100, nu = 4)) round(cbind(BPM = gamma.fit2$coef[sl.ind], Z = gamma.fit2$coef[sl.ind]/gamma.fit2$sd[sl.ind]), 4) ################################################################################# ## 5. Using interval time data ################################################################################# evlsint <- gen.evl(atus80int[, 1:3]) set.seed(570819) evlsint$eventlist <- evlsint$eventlist[sample(1:length(evlsint$eventlist), 500)] evlsint$event.key evlsint$eventlist$"20040706041558" int.base <- gen.intercepts(evlsint, type = 1, contr = FALSE) eT <- evlsint$event.key[, 2] supplist <- list() for (i in 1:length(evlsint$eventlist)) { supplist[[i]] <- matrix(0, nrow = NROW(evlsint$eventlist[[i]]), ncol = length(eT)) colnames(supplist[[i]]) <- eT } for (i in 1:length(evlsint$eventlist)) { evl <- evlsint$eventlist[[i]] supplist[[i]][1, grep("START", eT)] <- 1 for (j in 2:(NROW(evl))) { CE <- evl[j, 1] PE <- evl[j - 1, 1] if (PE == 0) { supplist[[i]][j, grep("STOP", eT)] <- 1 } if (PE != 0) { if (grepl("START", eT[PE])) { supplist[[i]][j, eT[CE]] <- 1 } if (grepl("STOP", eT[PE])) { supplist[[i]][j, grepl("START", eT)] <- 1 } } } } intfit.base <- rem(evlsint$eventlist, int.base, supplist = supplist, timing = "interval", estimator = "BPM") summary(intfit.base) bpm.start <- (exp(intfit.base$coef[grep("START", names(intfit.base$coef))])/sum(exp(intfit.base$coef[grep("START", names(intfit.base$coef))]))) in.ids <- which(atus80int$TUCASEID %in% names(evlsint$eventlist) & grepl("START", atus80int$Events)) mle.start <- prop.table(table(atus80int[in.ids, "Events"])) round(cbind(BPM = sort(bpm.start), MLE = sort(mle.start)), 4) sleepA <- paste("ab", c(letters[seq(1, 26, 2)], "A"), c(letters[seq(2, 26, 2)], "B"), "a", sep = "") sleepB <- paste("ab", c(letters[seq(1, 26, 2)], "A"), c(letters[seq(2, 26, 2)], "B"), "ab", sep = "") sleep.glbs <- glb.sformlist(evlsint, list(sleepA, sleepB, c("abefa", "abopa"), c("abefab", "abopab")), cond = TRUE, new.names = c("Inter.Sl.Ons", "Inter.Sl.Dur", "PerCare.Sl.Ons", "PerCare.Sl.Dur")) sleep.int <- slbind(sleep.glbs, int.base) intfit.sleep <- rem(evlsint$eventlist, sleep.int, supplist = supplist, timing = "interval", estimator = "BPM") summary(intfit.sleep) exp(-0.890707 + 4.989059)/sum(exp(intfit.sleep$coef[grep("START", names(intfit.sleep$coef))])) ################################################################################# ## A. Appendix ################################################################################# set.seed(867) n <- 5 tmax <- 25 roweff <- rnorm(n) roweff <- roweff - roweff[1] coleff <- rnorm(n) coleff <- coleff - coleff[1] lambda <- exp(outer(roweff, coleff, "+")) diag(lambda) <- 0 ratesum <- sum(lambda) esnd <- as.vector(row(lambda)) erec <- as.vector(col(lambda)) time <- 0 edgelist <- vector() while (time < tmax) { drawsr <- sample(1:(n^2), 1, prob = as.vector(lambda)) time <- time + rexp(1, ratesum) if (time < tmax) edgelist <- rbind(edgelist, c(time, esnd[drawsr], erec[drawsr])) else edgelist <- rbind(edgelist, c(tmax, NA, NA)) } evl <- gen.evl(cbind(apply(edgelist[, 2:3], 1, paste, collapse = "."), rep(1, NROW(edgelist))), null.events = "NA.NA") evl$event.key evl$event.key <- rbind(evl$event.key, c("u", "3.4")) ev.ints <- gen.intercepts(evl, contr = FALSE) sformlist <- c(lapply(2:n, function(z) { str <- paste(z, ".", sep = "") grep(str, evl$event.key[-20, 2]) }), lapply(2:n, function(z) { str <- paste(".", z, sep = "") grep(str, evl$event.key[-20, 2]) })) b1 <- lapply(sformlist, function(x) ev.ints[[1]][[1]][, , x]) b1.l <- lapply(b1, apply, MARGIN = 1:2, sum) FEs <- array(unlist(b1.l), dim = c(nrow(b1.l[[1]]), ncol(b1.l[[1]]), length(b1.l))) dimnames(FEs) <- list(dimnames(b1[[1]])[[1]], dimnames(b1[[1]])[[2]], c(paste("FESnd", 2:n, sep = "."), paste("FERec", 2:n, sep = "."))) fit.rem.dyad <- rem.dyad(edgelist, n, effects = c("FESnd", "FERec"), fit.method = "MLE", hessian = TRUE) fit.rem <- rem(evl$eventlist, sfl2statslist(list(FEs)), estimator = "MLE") lapply(list(rem.dyad = fit.rem.dyad, rem = fit.rem), summary) -2 * 813 * log(1/20) evl2 <- gen.evl(cbind(apply(edgelist[, 2:3], 1, paste, collapse = "."), edgelist[, 1], rep(1, NROW(edgelist))), null.events = "NA.NA") evl2$event.key <- rbind(evl2$event.key, c("u", "3.4")) fit.rem.dyad2 <- rem.dyad(edgelist, n, effects = c("FESnd", "FERec"), fit.method = "MLE", ordinal = FALSE, hessian = TRUE) fit.rem2 <- rem(evl2$eventlist, sfl2statslist(list(FEs)), estimator = "MLE", timing = "interval") lapply(list(rem.dyad = fit.rem.dyad2, rem = fit.rem2), summary)