Techniques of Event History Modeling, in R
- 2. Event History Data Structures
- 3. Descriptive Methods
- 4. Exponential Transition Rate Models
- 5. Piecewise Constant Exponential Models
- 6. Exponential Models with Time-dependent Covariates
- 7. Parametric Models of Time-dependence
- 8. Methods to Check Parametric Assumptions
- 9. Semi-parametric Transition Rate Models
- 10. Models with a Mixture Distribution
The command files distributed with Blossfeld and Rohwer,
Techniques of Event History Modeling, are the second body of
TDA examples, and the one most people met the program through. They live
in examples/ehhnew, about ninety of them, in eight
families:
| prefix | book chapter | what it covers |
|---|---|---|
ehb |
2 | episode data |
ehc |
3 | life tables and product-limit estimation |
ehd |
4 | exponential models |
ehe |
5 | piecewise constant models |
ehf |
6 | episode splitting and time-varying covariates |
ehg |
7 | parametric models |
ehh |
8 | graphical tests of distributional assumptions |
ehi |
9 | Cox models |
The suffixes are worth knowing, since they are not what they look
like: a file ending in a plots the fitted rate, while
s and p are the second and third steps of the
pseudoresidual check of the book’s section 8.2 – s
estimates a product-limit survivor function for a model’s residuals,
p plots log S against them. They are separate
command files in TDA only because each step writes its result to disk
for the next to read; here they follow the fit directly.
This vignette walks the families in order, naming the command file
each box corresponds to and, in brackets, the page and box number in the
book where the same output is printed, so the two can be read side by
side. tests/vignette-qa/vignette-qa-book.csv in the sources
is the other direction: one row per box and figure in the book, with the
vignette box that reproduces it. The data are the ones of the manual’s
section 3.3.3, and the covariate set is the book’s own.
Box 1The data, and the book’s covariates
d <- tda_rrdat() # rrdat.1, with the derived variables
d4 <- tda_rrdat(states = 4) # job changes split by prestige
XVARS <- c("EDU", "COHO2", "COHO3", "LFX", "PNOJ", "PRES")
with_x <- function(lhs) stats::reformulate(XVARS, response = lhs)Most plots below share one coordinate system, taken from the command
files’ own psetup rather than from the range of the data,
so that a picture can be laid beside the book.
2. Event History Data Structures
Book chapter 2, pp. 41-55. Command files ehb1 to
ehb3.
ehb1 is the data before any episode definition, and ehb2 the episode data itself: one transition, from 0 to 1.
Box 2ehb1, ehb2: the data and the episodes [book pp. 43-52, Boxes 2.2.1-2.2.5]
# the twelve variables of rrdat.1; DES and TFP are derived by ehb1.cf
str(d[, 1:12])'data.frame': 600 obs. of 12 variables:
$ ID : int 1 2 2 2 3 3 3 3 3 4 ...
$ NOJ : int 1 1 2 3 1 2 3 4 5 1 ...
$ TStart: int 555 593 639 673 688 700 730 742 817 872 ...
$ TFin : int 982 638 672 892 699 729 741 816 828 926 ...
$ SEX : int 1 2 2 2 2 2 2 2 2 2 ...
$ TI : int 982 982 982 982 982 982 982 982 982 982 ...
$ TB : int 351 357 357 357 473 473 473 473 473 604 ...
$ TE : int 555 593 593 593 688 688 688 688 688 872 ...
$ TMAR : int 679 762 762 762 870 870 870 870 870 872 ...
$ PRES : int 34 22 46 46 41 41 44 44 44 55 ...
$ PRESN : int -1 46 46 -1 41 44 44 44 -1 -1 ...
$ EDU : int 17 10 10 10 11 11 11 11 11 13 ...head(d[, 1:12], 9) ID NOJ TStart TFin SEX TI TB TE TMAR PRES PRESN EDU
1 1 1 555 982 1 982 351 555 679 34 -1 17
2 2 1 593 638 2 982 357 593 762 22 46 10
3 2 2 639 672 2 982 357 593 762 46 46 10
4 2 3 673 892 2 982 357 593 762 46 -1 10
5 3 1 688 699 2 982 473 688 870 41 41 11
6 3 2 700 729 2 982 473 688 870 41 44 11
7 3 3 730 741 2 982 473 688 870 44 44 11
8 3 4 742 816 2 982 473 688 870 44 44 11
9 3 5 817 828 2 982 473 688 870 44 -1 11head(d[, c("ID", "NOJ", "TStart", "TFin", "SEX", "TFP", "DES")], 4) ID NOJ TStart TFin SEX TFP DES
1 1 1 555 982 1 428 0
2 2 1 593 638 2 46 1
3 2 2 639 672 2 34 1
4 2 3 673 892 2 220 1tda_rate(Surv(TFP, DES) ~ 1, d)$episodes sn org des episodes weighted mean.duration ts.min tf.max excluded
1 1 0 0 142 142 128.17606 0 428 FALSE
2 1 0 1 458 458 49.30349 0 350 FALSEehb3 is the multi-episode form of the same data (Box
2.2.6): the spells of a person are chained, the origin state of spell k
is 1 + 2(k - 1), the destination the next number up if the
job ended, and time runs from the person’s first job. id=
and spell= are edef(id=, sn=); the episode
table is the book’s Box 2.2.8.
Box 3ehb3: multi-episode data [book pp. 53-55, Boxes 2.2.6-2.2.8]
dm <- d
dm$ORG <- 1 + 2 * (dm$NOJ - 1)
dm$DESM <- ifelse(dm$TFin < dm$TI, dm$ORG + 1, dm$ORG)
dm$TFC <- ave(dm$TStart, dm$ID, FUN = function(x) x[1])
dm$TSP <- dm$TStart - dm$TFC
dm$TFPM <- dm$TFin - dm$TFC + 1
# the episodes of the first three individuals, as edef() will see them
head(dm[dm$ID <= 3, c("ID", "NOJ", "TStart", "TFin", "SEX", "ORG", "DESM",
"TSP", "TFPM")], 9) ID NOJ TStart TFin SEX ORG DESM TSP TFPM
1 1 1 555 982 1 1 1 0 428
2 2 1 593 638 2 1 2 0 46
3 2 2 639 672 2 3 4 46 80
4 2 3 673 892 2 5 6 80 300
5 3 1 688 699 2 1 2 0 12
6 3 2 700 729 2 3 4 12 42
7 3 3 730 741 2 5 6 42 54
8 3 4 742 816 2 7 8 54 129
9 3 5 817 828 2 9 10 129 141eb3 <- suppressWarnings(
tda_rate(Surv(TSP, TFPM, ORG, DESM) ~ 1, dm, id = "ID", spell = "NOJ"))
eb3$episodes sn org des episodes weighted mean.duration ts.min tf.max excluded
1 1 1 1 16 16 233.43750 0 428 FALSE
2 1 1 2 185 185 52.72973 0 326 FALSE
3 2 3 3 36 36 150.16667 13 450 FALSE
4 2 3 4 126 126 52.30159 3 406 FALSE
5 3 5 5 38 38 121.73684 9 463 FALSE
6 3 5 6 69 69 50.36232 8 403 FALSE
7 4 7 7 24 24 85.33334 52 465 FALSE
8 4 7 8 38 38 37.71053 9 329 FALSE
9 5 9 9 9 9 110.22222 66 441 FALSE
10 5 9 10 23 23 22.60870 39 348 FALSE
11 6 11 11 8 8 45.75000 151 379 FALSE
12 6 11 12 12 12 43.75000 46 324 FALSE
13 7 13 13 7 7 114.14286 132 427 FALSE
14 7 13 14 4 4 52.75000 112 368 FALSE
15 8 15 15 3 3 65.00000 129 440 FALSE
16 8 15 16 1 1 72.00000 213 285 FALSE
17 9 17 17 1 1 34.00000 285 319 TRUEattr(eb3$episodes, "sums") # the totals under each spell; sn = NA is all spells sn episodes weighted
1 1 201 201
2 2 162 162
3 3 107 107
4 4 62 62
5 5 32 32
6 6 20 20
7 7 11 11
8 8 4 4
9 9 1 1
10 NA 600 600The ninth spell holds a single, censored episode – no transition to
estimate – and TDA marks it with * in the Excl
column, which is excluded = TRUE here; the warning that
says so is the one suppressed above.
3. Descriptive Methods
Book chapter 3, pp. 57-85. Life tables and the product-limit
estimator; command files ehc1 to
ehc10.
ehc1 is the life table over 30-month intervals, ehc3 the same by sex.
Box 4ehc1, ehc3: life tables [book pp. 60-66, Boxes 3.1.1-3.1.4]
lt <- tda_ltb(Surv(TFP, DES) ~ 1, d, tp = seq(0, 500, by = 30))
summary(lt, n = nrow(lt$tables[[1]]))Call: tda_ltb(formula = Surv(TFP, DES) ~ 1, data = d, tp = seq(0, 500,
by = 30))
Life table. SN 1. Origin state 0.
start midpoint entering censored exposed events prob
0 15 600 28 586.0 223 0.38054608
30 45 349 23 337.5 113 0.33481481
60 75 213 15 205.5 51 0.24817518
90 105 147 16 139.0 25 0.17985612
120 135 106 15 98.5 24 0.24365482
150 165 67 5 64.5 9 0.13953488
180 195 53 9 48.5 4 0.08247423
210 225 40 5 37.5 3 0.08000000
240 255 32 5 29.5 0 0.00000000
270 285 27 7 23.5 2 0.08510638
300 315 18 5 15.5 2 0.12903226
330 345 11 1 10.5 2 0.19047619
360 375 8 3 6.5 0 0.00000000
390 405 5 4 3.0 0 0.00000000
420 435 1 1 0.5 0 0.00000000
Median duration: 47.27858
Cases: 600 weighted: 600lt2 <- tda_ltb(Surv(TFP, DES) ~ as.factor(SEX), d, tp = seq(0, 500, by = 30))
lt2$episodes # the episodes behind both groups sn org des episodes weighted mean.duration ts.min tf.max excluded
1 1 0 0 142 142 128.17606 0 428 FALSE
2 1 0 1 458 458 49.30349 0 350 FALSElt2$summary sn origin cases weighted median group
1 1 0 348 348 56.84110 1
2 1 0 252 252 39.63376 2ehc2 plots the survivor function the table implies, and ehc2a and ehc2b the density and the rate, each on its own scale.
Box 5ehc2: the life table survivor function [book p. 65, Fig. 3.1.1]
sv <- tda_survivor(lt)
p <- tda_ps(sv, width = 80, height = 40, xlim = c(0, 300), ylim = c(0, 1))
p <- tda_pl(p, "plxa", sc = 60, ic = 5)
p <- tda_pl(p, "plya", sc = 1, ic = 10)
p <- tda_pl_lines(tda_pl(p, "plframe"), "time", "survivor")
plot(tda_pl_text(p, "Life Table Survivor Function", at = c(80, 0.8)))
Box 6ehc2a, ehc2b: the density and the rate [book p. 65, Fig. 3.1.1]
p <- tda_ps(sv, width = 80, height = 40, xlim = c(0, 300), ylim = c(0, 0.02))
p <- tda_pl(p, "plxa", sc = 60, ic = 5)
p <- tda_pl(p, "plya", sc = 0.01, ic = 10)
p <- tda_pl_lines(tda_pl(p, "plframe"), "time", "density")
plot(tda_pl_text(p, "Life Table Density Function", at = c(80, 0.017)))
p <- tda_ps(sv, width = 80, height = 40, xlim = c(0, 300), ylim = c(0, 0.02))
p <- tda_pl(p, "plxa", sc = 60, ic = 5)
p <- tda_pl(p, "plya", sc = 0.01, ic = 10)
p <- tda_pl_lines(tda_pl(p, "plframe"), "time", "rate")
plot(tda_pl_text(p, "Life Table Transition Rate", at = c(80, 0.017)))
ehc4 draws one curve per sex.
Box 7ehc4: survivor functions by sex [book p. 67, Fig. 3.1.2]
sv2 <- tda_survivor(lt2)
p <- tda_ps(sv2, width = 80, height = 40, xlim = c(0, 300), ylim = c(0, 1))
p <- tda_pl(p, "plxa", sc = 60, ic = 5)
p <- tda_pl(p, "plya", sc = 1, ic = 10)
p <- tda_pl_lines(p, "time", "survivor", by = "group")
p <- tda_pl_text(p, "Men", at = c(240, 0.25))
p <- tda_pl_text(p, "Women", at = c(240, 0.09))
plot(tda_pl_text(p, "Life Table Survivor Function", at = c(80, 0.8)))
ehc5 is the product-limit estimate, ehc7 the comparison of the two curves, and ehc9 the four destination states.
Box 8ehc5, ehc7, ehc9: product-limit estimation [book pp. 74-83, Boxes 3.2.1-3.2.2, 3.3.1-3.3.4a, 3.3.4b]
pl <- tda_ple(Surv(TFP, DES) ~ 1, d)
pl$episodes # the episode table TDA prints above the estimate sn org des episodes weighted mean.duration ts.min tf.max excluded
1 1 0 0 142 142 128.17606 0 428 FALSE
2 1 0 1 458 458 49.30349 0 350 FALSEprint(pl, n = 10)Call: tda_ple(formula = Surv(TFP, DES) ~ 1, data = d)
SN 1. Transition: 0,1 - Product-Limit Estimation
id index time events censored n.risk survivor std.err cum.rate
0 0 0 0 0 600 1.0000000 0.000000000 0.000000000
0 1 2 2 0 600 0.9966667 0.002353091 0.003338901
0 2 3 5 1 597 0.9883194 0.004389012 0.011749380
0 3 4 9 2 590 0.9732433 0.006599178 0.027121160
0 4 5 3 0 581 0.9682180 0.007174619 0.032298048
0 5 6 10 1 577 0.9514378 0.008796359 0.049781011
0 6 7 9 0 567 0.9363356 0.009993927 0.065781352
0 7 8 6 1 557 0.9262494 0.010700980 0.076611783
0 8 9 7 3 548 0.9144177 0.011460673 0.089467791
0 9 10 8 1 540 0.9008708 0.012250884 0.104393441
... 122 more rows
Median duration: 43.03101
Duration times limited to: 350
Cases: 600 weighted: 600pl2 <- tda_ple(Surv(TFP, DES) ~ as.factor(SEX), d, compare = TRUE)
pl2$comparison sn org des test statistic df signif
1 1 0 1 Log-Rank (Savage) 20.603209 1 0.9999943
2 1 0 1 Wilcoxon (Breslow) 9.364723 1 0.9977880
3 1 0 1 Wilcoxon (Tarone-Ware) 14.326703 1 0.9998463
4 1 0 1 Wilcoxon (Prentice) 10.695087 1 0.9989258# one product-limit table per destination: upward, then lateral moves
pl4 <- tda_ple(Surv(TFP, DES) ~ 1, d4)
pl4$episodes sn org des episodes weighted mean.duration ts.min tf.max excluded
1 1 0 0 142 142 128.17606 0 428 FALSE
2 1 0 1 84 84 45.04762 0 326 FALSE
3 1 0 2 219 219 45.91324 0 350 FALSE
4 1 0 3 155 155 56.40000 0 332 FALSEpl4$summary transition median limit cases weighted
1 0,1 NA NA 600 600
2 0,2 122.2409 350 600 600
3 0,3 215.9553 332 600 600for (b in pl4$blocks[1:2]) print(utils::head(b, 10), row.names = FALSE) id index time events censored n.risk survivor std.err cum.rate
0 0 0 0 0 600 1.0000000 0.000000000 0.000000000
0 1 2 1 0 600 0.9983333 0.001665277 0.001668057
0 2 4 2 10 590 0.9949492 0.002908833 0.005063646
0 3 6 3 10 577 0.9897761 0.004152988 0.010276516
0 4 7 1 9 567 0.9880305 0.004497599 0.012041742
0 5 8 1 11 557 0.9862566 0.004826665 0.013838687
0 6 9 1 7 548 0.9844569 0.005142460 0.015665172
0 7 10 3 6 540 0.9789877 0.006005599 0.021236217
0 8 12 5 15 524 0.9696462 0.007257296 0.030824018
0 9 13 1 18 499 0.9677030 0.007498389 0.032830037
id index time events censored n.risk survivor std.err cum.rate
1 0 0 0 0 600 1.0000000 0.000000000 0.000000000
1 1 3 2 3 597 0.9966499 0.002364896 0.003355708
1 2 4 4 5 590 0.9898930 0.004105345 0.010158455
1 3 5 2 5 581 0.9864854 0.004745918 0.013606734
1 4 6 4 5 577 0.9796467 0.005815798 0.020563284
1 5 7 3 3 567 0.9744634 0.006509588 0.025868336
1 6 8 4 7 557 0.9674654 0.007343240 0.033075575
1 7 9 5 5 548 0.9586382 0.008269553 0.042241542
1 8 10 2 4 540 0.9550877 0.008611598 0.045952121
1 9 11 3 9 528 0.9496611 0.009114803 0.051650142ehc6 plots the estimate, ehc10 one curve per destination state, and ehc8 adds the confidence bands, cut at 290 months where they widen past usefulness.
Box 9ehc6, ehc10: the estimated survivor functions [book pp. 76-85, Box 3.2.3; Fig. 3.2.1, Fig. 3.2.3]
p <- tda_ps(tda_survivor(pl), width = 80, height = 40, xlim = c(0, 300), ylim = c(0, 1))
p <- tda_pl(p, "plxa", sc = 60, ic = 5)
p <- tda_pl(p, "plya", sc = 1, ic = 10)
p <- tda_pl_lines(p, "time", "survivor")
plot(tda_pl_text(p, "Product-Limit Survivor Function", at = c(60, 0.8)))
s4 <- tda_survivor(pl4)
p <- tda_ps(s4, width = 80, height = 40, xlim = c(0, 300), ylim = c(0, 1))
p <- tda_pl(p, "plxa", sc = 60, ic = 5)
p <- tda_pl(p, "plya", sc = 1, ic = 10)
for (i in seq_along(g4 <- c("0,1", "0,2", "0,3")))
p <- tda_pl_lines(p, "time", "survivor", by = "group",
select = g4[i], lty = c(1, 5, 8)[i])
p <- tda_pl_text(p, "upward moves", at = c(200, 0.75))
p <- tda_pl_text(p, "lateral moves", at = c(200, 0.3))
plot(tda_pl_text(p, "downward moves", at = c(200, 0.55)))
A band is two fills, the upper curve to the axis in grey and the
lower curve to the axis in white on top of it; ehc8.cf
draws the axes first, so the white fill lies over half the width of the
x axis line and the axis looks thinner under the bands than beyond them.
The axes go last here.
Box 10ehc8: with confidence bands [book p. 78, Fig. 3.3.1]
sb <- tda_survivor(pl2)
sb <- sb[sb$time < 290, ]
p <- tda_ps(sb, width = 80, height = 40, xlim = c(0, 300), ylim = c(0, 1))
p <- tda_pl_lines(p, "time", "survivor", by = "group", select = "1",
band = c("upper", "lower"), lty = 1)
p <- tda_pl_lines(p, "time", "survivor", by = "group", select = "2",
band = c("upper", "lower"), lty = 5)
p <- tda_pl(p, "plxa", sc = 60, ic = 5)
p <- tda_pl(p, "plya", sc = 1, ic = 10)
p <- tda_pl_text(p, "Men", at = c(240, 0.30))
p <- tda_pl_text(p, "Women", at = c(50, 0.10))
plot(tda_pl_text(p, "Survivor Functions for Men and Women", at = c(60, 0.8)))
4. Exponential Transition Rate Models
Book chapter 4, pp. 87-118. Command files ehd1 to
ehd7.
ehd1 is the null model and ehd2 the same with the book’s covariates.
Box 11ehd1, ehd2: exponential models [book pp. 91-98, Boxes 4.1.1-4.1.5]
ehd1 <- tda_rate(Surv(TFP, DES) ~ 1, d, prate = "0(10)300")
ehd1Call: tda_rate(formula = Surv(TFP, DES) ~ 1, data = d, prate = "0(10)300")
Episodes: 600
Model: exponential (TDA 2)
logLik (starting values): -2514.02
logLik: -2514.02
Converged in 1 iterations
Idx SN Org Des MT Variable Coeff Error C/Error Signif
1 1 0 1 A Constant -4.4891 0.0467 -96.0715 1.0000ehd1$episodes # the episodes TDA built, before the estimates sn org des episodes weighted mean.duration ts.min tf.max excluded
1 1 0 0 142 142 128.17606 0 428 FALSE
2 1 0 1 458 458 49.30349 0 350 FALSEtda_rates(ehd1) group time rate survivor density
1 0 0 0.01123044 1.00000000 0.0112304448
2 0 10 0.01123044 0.89377211 0.0100374583
3 0 20 0.01123044 0.79882858 0.0089712003
4 0 30 0.01123044 0.71397071 0.0080182086
5 0 40 0.01123044 0.63812710 0.0071664512
6 0 50 0.01123044 0.57034021 0.0064051742
7 0 60 0.01123044 0.50975417 0.0057247661
8 0 70 0.01123044 0.45560406 0.0051166362
9 0 80 0.01123044 0.40720620 0.0045731068
10 0 90 0.01123044 0.36394955 0.0040873153
11 0 100 0.01123044 0.32528795 0.0036531284
12 0 110 0.01123044 0.29073330 0.0032650643
13 0 120 0.01123044 0.25984931 0.0029182234
14 0 130 0.01123044 0.23224607 0.0026082267
15 0 140 0.01123044 0.20757506 0.0023311602
16 0 150 0.01123044 0.18552480 0.0020835260
17 0 160 0.01123044 0.16581689 0.0018621974
18 0 170 0.01123044 0.14820251 0.0016643801
19 0 180 0.01123044 0.13245927 0.0014875765
20 0 190 0.01123044 0.11838840 0.0013295544
21 0 200 0.01123044 0.10581225 0.0011883187
22 0 210 0.01123044 0.09457204 0.0010620861
23 0 220 0.01123044 0.08452585 0.0009492629
24 0 230 0.01123044 0.07554685 0.0008484247
25 0 240 0.01123044 0.06752167 0.0007582983
26 0 250 0.01123044 0.06034898 0.0006777459
27 0 260 0.01123044 0.05393824 0.0006057504
28 0 270 0.01123044 0.04820849 0.0005414028
29 0 280 0.01123044 0.04308741 0.0004838907
30 0 290 0.01123044 0.03851032 0.0004324880
31 0 300 0.01123044 0.03441945 0.0003865457summary(tda_rate(with_x("Surv(TFP, DES)"), d))Call: tda_rate(formula = with_x("Surv(TFP, DES)"), data = d)
Episodes: 600
Model: exponential
Converged in 5 iterations
Estimate Std. Error C/Error Pr(>|t|)
Constant -4.48944489 0.27950029 -16.0624 < 2.2e-16 ***
EDU 0.07730132 0.02470334 3.1292 0.0017529 **
COHO2 0.60803578 0.11355087 5.3547 8.568e-08 ***
COHO3 0.61079970 0.11854198 5.1526 2.569e-07 ***
LFX -0.00317926 0.00093783 -3.3900 0.0006989 ***
PNOJ 0.05963838 0.04415253 1.3507 0.1767802
PRES -0.02800646 0.00553034 -5.0641 4.102e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
logLik (starting values): -2514.02
logLik: -2465.987 df: 7 AIC: 4945.975
Null model (exponential, constant rate) logLik: -2514.02
LR test against it: chi2 = 96.0657 on 6 df, p = < 2.2e-16The rate this model implies is a constant, exp(-4.4891)
= 0.0112, and the book prints it as the table ehd1.prs
above rather than drawing it. Everything that follows is an attempt to
let that constant vary.
ehd3 estimates the same specification over three destination states.
Box 12ehd3: three destinations [book pp. 102-104, Boxes 4.2.1-4.2.2]
ehd3 <- tda_rate(with_x("Surv(TFP, DES)"), d4)
ehd3$episodes sn org des episodes weighted mean.duration ts.min tf.max excluded
1 1 0 0 142 142 128.17606 0 428 FALSE
2 1 0 1 84 84 45.04762 0 326 FALSE
3 1 0 2 219 219 45.91324 0 350 FALSE
4 1 0 3 155 155 56.40000 0 332 FALSEehd3Call: tda_rate(formula = with_x("Surv(TFP, DES)"), data = d4)
Episodes: 600
Model: exponential (TDA 2)
logLik (starting values): -2986
logLik: -2884.788
Converged in 7 iterations
Transitions: 0->1, 0->2, 0->3
Idx SN Org Des MT Variable Coeff Error C/Error Signif
1 1 0 1 A Constant -5.1036 0.6175 -8.2648 1.0000
2 1 0 1 A EDU 0.3018 0.0430 7.0215 1.0000
3 1 0 1 A COHO2 0.6379 0.2713 2.3510 0.9813
4 1 0 1 A COHO3 0.7354 0.2765 2.6595 0.9922
5 1 0 1 A LFX -0.0023 0.0021 -1.0898 0.7242
6 1 0 1 A PNOJ 0.1731 0.1004 1.7239 0.9153
7 1 0 1 A PRES -0.1441 0.0141 -10.1970 1.0000
8 1 0 2 A Constant -5.8043 0.4054 -14.3164 1.0000
9 1 0 2 A EDU 0.0033 0.0380 0.0881 0.0702
10 1 0 2 A COHO2 0.6722 0.1643 4.0926 1.0000
11 1 0 2 A COHO3 0.6843 0.1732 3.9501 0.9999
12 1 0 2 A LFX -0.0031 0.0014 -2.2579 0.9760
13 1 0 2 A PNOJ 0.0323 0.0645 0.5010 0.3836
14 1 0 2 A PRES 0.0081 0.0081 1.0094 0.6872
15 1 0 3 A Constant -4.9051 0.5251 -9.3405 1.0000
16 1 0 3 A EDU -0.0167 0.0490 -0.3404 0.2664
17 1 0 3 A COHO2 0.5300 0.1935 2.7387 0.9938
18 1 0 3 A COHO3 0.5705 0.2030 2.8101 0.9950
19 1 0 3 A LFX -0.0037 0.0017 -2.2649 0.9765
20 1 0 3 A PNOJ 0.0407 0.0775 0.5252 0.4005
21 1 0 3 A PRES -0.0157 0.0097 -1.6208 0.8949ehd4 and ehd5 take random
subsamples so that the destinations are comparably sized.
tda_runif() is TDA’s generator, a multiplicative
congruential one with a fixed seed, so the episodes drawn here are the
episodes TDA drew; runif() would give a different subsample
and so different estimates.
Box 13ehd4, ehd5: subsamples, one transition at a time [book pp. 106-110, Boxes 4.2.3-4.2.6]
sub <- d4[tda_runif(nrow(d4)) <= 84 / 219, ]
sub$DES <- ifelse(sub$DES == 2, 2L, 0L)
cat(sprintf("%d of %d episodes selected\n", nrow(sub), nrow(d4)))235 of 600 episodes selectedehd4 <- tda_rate(with_x("Surv(TFP, DES)"), sub)
ehd4$episodes sn org des episodes weighted mean.duration ts.min tf.max excluded
1 1 0 0 145 145 82.71034 0 428 FALSE
2 1 0 2 90 90 49.33333 0 209 FALSEehd4Call: tda_rate(formula = with_x("Surv(TFP, DES)"), data = sub)
Episodes: 235
Model: exponential (TDA 2)
logLik (starting values): -558.6513
logLik: -552.8855
Converged in 5 iterations
Idx SN Org Des MT Variable Coeff Error C/Error Signif
1 1 0 2 A Constant -5.3301 0.6960 -7.6583 1.0000
2 1 0 2 A EDU -0.0097 0.0594 -0.1627 0.1293
3 1 0 2 A COHO2 0.6587 0.2516 2.6181 0.9912
4 1 0 2 A COHO3 0.4869 0.2732 1.7824 0.9253
5 1 0 2 A LFX -0.0025 0.0023 -1.1042 0.7305
6 1 0 2 A PNOJ -0.0109 0.1138 -0.0961 0.0765
7 1 0 2 A PRES 0.0014 0.0127 0.1073 0.0854d5 <- d4[tda_runif(nrow(d4)) <= 84 / 155, ]
d5$DES <- ifelse(d5$DES == 3, 3L, 0L)
cat(sprintf("%d of %d episodes selected\n", nrow(d5), nrow(d4)))314 of 600 episodes selectedehd5 <- tda_rate(with_x("Surv(TFP, DES)"), d5)
ehd5$episodes sn org des episodes weighted mean.duration ts.min tf.max excluded
1 1 0 0 235 235 76.25106 0 428 FALSE
2 1 0 3 79 79 50.83544 0 312 FALSEehd5Call: tda_rate(formula = with_x("Surv(TFP, DES)"), data = d5)
Episodes: 314
Model: exponential (TDA 2)
logLik (starting values): -523.4849
logLik: -516.5031
Converged in 6 iterations
Idx SN Org Des MT Variable Coeff Error C/Error Signif
1 1 0 3 A Constant -5.6720 0.7092 -7.9981 1.0000
2 1 0 3 A EDU 0.0169 0.0601 0.2806 0.2210
3 1 0 3 A COHO2 0.5486 0.2714 2.0218 0.9568
4 1 0 3 A COHO3 0.5347 0.2890 1.8498 0.9357
5 1 0 3 A LFX -0.0056 0.0027 -2.0972 0.9640
6 1 0 3 A PNOJ 0.1471 0.1122 1.3114 0.8103
7 1 0 3 A PRES -0.0085 0.0131 -0.6488 0.4836Box 4.2.7 of the book lays the two sets of estimates side by side,
which is the point of the exercise: with the events standardized,
COHO3 and LFX lose their significance for
lateral moves and COHO3 for downward ones, so the same data
can support different substantive conclusions depending on how many
events each transition happens to have. Upward moves are the baseline
and appear in both columns unchanged.
Box 14Box 4.2.7: standardized against unstandardized estimates [book p. 111, Box 4.2.7]
unstd <- function(des) {
e <- ehd3$estimates
e[e$Des == des, c("Variable", "Coeff", "C/Error")]
}
std <- list(`1` = unstd(1),
`2` = ehd4$estimates[, c("Variable", "Coeff", "C/Error")],
`3` = ehd5$estimates[, c("Variable", "Coeff", "C/Error")])
for (des in 1:3) {
a <- unstd(des)
b <- std[[as.character(des)]]
cmp <- data.frame(Variable = a$Variable,
Coeff = a$Coeff, `C/Error` = a$`C/Error`,
Coeff.s = b$Coeff, `C/Error.s` = b$`C/Error`,
check.names = FALSE)
cat(sprintf("\n%s moves. Events: %d (unstandardized) %d (standardized)\n",
c("Upward", "Lateral", "Downward")[des],
ehd3$episodes$episodes[ehd3$episodes$des == des],
if (des == 1) ehd3$episodes$episodes[ehd3$episodes$des == 1]
else list(NULL, ehd4, ehd5)[[des]]$episodes$episodes[2]))
print(format(cmp, digits = 4), row.names = FALSE)
}
Upward moves. Events: 84 (unstandardized) 84 (standardized)
Variable Coeff C/Error Coeff.s C/Error.s
Constant -5.103617 -8.265 -5.103617 -8.265
EDU 0.301790 7.022 0.301790 7.022
COHO2 0.637895 2.351 0.637895 2.351
COHO3 0.735415 2.660 0.735415 2.660
LFX -0.002265 -1.090 -0.002265 -1.090
PNOJ 0.173106 1.724 0.173106 1.724
PRES -0.144116 -10.197 -0.144116 -10.197
Lateral moves. Events: 219 (unstandardized) 90 (standardized)
Variable Coeff C/Error Coeff.s C/Error.s
Constant -5.804317 -14.31637 -5.330081 -7.65834
EDU 0.003345 0.08806 -0.009672 -0.16271
COHO2 0.672225 4.09264 0.658747 2.61809
COHO3 0.684335 3.95009 0.486906 1.78241
LFX -0.003085 -2.25785 -0.002496 -1.10420
PNOJ 0.032293 0.50100 -0.010935 -0.09608
PRES 0.008136 1.00939 0.001362 0.10729
Downward moves. Events: 155 (unstandardized) 79 (standardized)
Variable Coeff C/Error Coeff.s C/Error.s
Constant -4.905088 -9.3405 -5.671986 -7.9981
EDU -0.016686 -0.3404 0.016866 0.2806
COHO2 0.529999 2.7387 0.548629 2.0218
COHO3 0.570488 2.8101 0.534684 1.8498
LFX -0.003745 -2.2649 -0.005638 -2.0972
PNOJ 0.040694 0.5252 0.147083 1.3114
PRES -0.015691 -1.6208 -0.008525 -0.6488ehd6 holds the cohort effects equal across the three transitions, which is what a constraint is for.
Box 15ehd6: cohort effects constrained equal [book pp. 112-113, Boxes 4.2.8-4.2.9]
ehd6 <- tda_rate(with_x("Surv(TFP, DES)"), d4,
constraints = c("b3 - b10 = 0", "b3 - b17 = 0",
"b4 - b11 = 0", "b4 - b18 = 0"))
ehd6Call: tda_rate(formula = with_x("Surv(TFP, DES)"), data = d4, constraints = c("b3 - b10 = 0",
"b3 - b17 = 0", "b4 - b11 = 0", "b4 - b18 = 0"))
Episodes: 600
Model: exponential (TDA 2)
logLik (starting values): -2986
logLik: -2885.017
Converged in 7 iterations
Transitions: 0->1, 0->2, 0->3
Idx SN Org Des MT Variable Coeff Error C/Error Signif
1 1 0 1 A Constant -5.0782 0.5981 -8.4902 1.0000
2 1 0 1 A EDU 0.3006 0.0423 7.1066 1.0000
3 1 0 1 A COHO2 0.6182 0.1136 5.4404 1.0000
4 1 0 1 A COHO3 0.6543 0.1190 5.4985 1.0000
5 1 0 1 A LFX -0.0024 0.0020 -1.1832 0.7633
6 1 0 1 A PNOJ 0.1780 0.0991 1.7962 0.9275
7 1 0 1 A PRES -0.1435 0.0140 -10.2516 1.0000
8 1 0 2 A Constant -5.7656 0.3950 -14.5984 1.0000
9 1 0 2 A EDU 0.0017 0.0376 0.0445 0.0355
10 1 0 2 A COHO2 0.6182 0.1136 5.4404 1.0000
11 1 0 2 A COHO3 0.6543 0.1190 5.4985 1.0000
12 1 0 2 A LFX -0.0031 0.0013 -2.3396 0.9807
13 1 0 2 A PNOJ 0.0341 0.0641 0.5328 0.4058
14 1 0 2 A PRES 0.0083 0.0080 1.0365 0.7000
15 1 0 3 A Constant -4.9621 0.5197 -9.5479 1.0000
16 1 0 3 A EDU -0.0152 0.0491 -0.3092 0.2428
17 1 0 3 A COHO2 0.6182 0.1136 5.4404 1.0000
18 1 0 3 A COHO3 0.6543 0.1190 5.4985 1.0000
19 1 0 3 A LFX -0.0036 0.0016 -2.2089 0.9728
20 1 0 3 A PNOJ 0.0362 0.0774 0.4683 0.3604
21 1 0 3 A PRES -0.0161 0.0097 -1.6623 0.9035ehd6$episodes # the episodes TDA built, as it reports them sn org des episodes weighted mean.duration ts.min tf.max excluded
1 1 0 0 142 142 128.17606 0 428 FALSE
2 1 0 1 84 84 45.04762 0 326 FALSE
3 1 0 2 219 219 45.91324 0 350 FALSE
4 1 0 3 155 155 56.40000 0 332 FALSEehd7 keeps the first four spells and gives each its
coefficients: the effect of education on leaving a first job need not be
its effect on leaving a fourth. id= and spell=
are TDA’s edef(id=, sn=), and spell_terms=
gives one spell a covariate set of its own – the first job has no labour
force experience to speak of, so LFX is left out of it, as
ehd7.cf does.
Box 16ehd7: a multi-episode model [book pp. 115-118, Boxes 4.3.1-4.3.2a, 4.3.2b]
d7 <- d[d$NOJ <= 4, ]
d7 <- d7[order(d7$ID, d7$NOJ), ]
# ehd7.cf measures time from the start of the FIRST job, not from the
# start of each spell: TFC is the first spell's starting time
# carried forward, and the episode runs from TStart - TFC to
# TFin - TFC + 1. Using each spell's duration instead would start
# every spell at 0 and give a different episode structure.
d7$TFC <- ave(ifelse(d7$NOJ == 1, d7$TStart, NA), d7$ID,
FUN = function(v) v[!is.na(v)][1])
d7$TSP <- d7$TStart - d7$TFC
d7$TFPC <- d7$TFin - d7$TFC + 1
ehd7 <- tda_rate(Surv(TSP, TFPC, DES) ~ EDU + COHO2 + COHO3 + LFX + PRES, d7,
id = "ID", spell = "NOJ",
spell_terms = list("1" = c("EDU", "COHO2", "COHO3", "PRES")))
ehd7$episodes sn org des episodes weighted mean.duration ts.min tf.max excluded
1 1 0 0 16 16 233.43750 0 428 FALSE
2 1 0 1 185 185 52.72973 0 326 FALSE
3 2 0 0 36 36 150.16667 13 450 FALSE
4 2 0 1 126 126 52.30159 3 406 FALSE
5 3 0 0 38 38 121.73684 9 463 FALSE
6 3 0 1 69 69 50.36232 8 403 FALSE
7 4 0 0 24 24 85.33334 52 465 FALSE
8 4 0 1 38 38 37.71053 9 329 FALSEehd7Call: tda_rate(formula = Surv(TSP, TFPC, DES) ~ EDU + COHO2 + COHO3 +
LFX + PRES, data = d7, id = "ID", spell = "NOJ", spell_terms = list(`1` = c("EDU",
"COHO2", "COHO3", "PRES")))
Episodes: 532
Model: exponential (TDA 2)
logLik (starting values): -2286.085
logLik: -2233.703
Converged in 6 iterations
Idx SN Org Des MT Variable Coeff Error C/Error Signif
1 1 0 1 A Constant -5.0179 0.4174 -12.0225 1.0000
2 1 0 1 A EDU 0.0622 0.0327 1.9052 0.9432
3 1 0 1 A COHO2 0.7221 0.1888 3.8245 0.9999
4 1 0 1 A COHO3 0.6638 0.1792 3.7050 0.9998
5 1 0 1 A PRES -0.0095 0.0091 -1.0452 0.7041
6 2 0 1 A Constant -3.7352 0.5778 -6.4644 1.0000
7 2 0 1 A EDU 0.1261 0.0516 2.4442 0.9855
8 2 0 1 A COHO2 0.3423 0.2157 1.5867 0.8874
9 2 0 1 A COHO3 0.2704 0.2334 1.1582 0.7532
10 2 0 1 A LFX -0.0062 0.0018 -3.3525 0.9992
11 2 0 1 A PRES -0.0510 0.0110 -4.6237 1.0000
12 3 0 1 A Constant -4.2628 0.7397 -5.7631 1.0000
13 3 0 1 A EDU 0.0522 0.0749 0.6976 0.5146
14 3 0 1 A COHO2 0.1975 0.2972 0.6645 0.4936
15 3 0 1 A COHO3 0.7143 0.3289 2.1719 0.9701
16 3 0 1 A LFX -0.0039 0.0020 -1.9843 0.9528
17 3 0 1 A PRES -0.0223 0.0139 -1.6008 0.8906
18 4 0 1 A Constant -2.2232 1.3489 -1.6481 0.9007
19 4 0 1 A EDU 0.0438 0.1329 0.3295 0.2583
20 4 0 1 A COHO2 1.1354 0.3830 2.9643 0.9970
21 4 0 1 A COHO3 0.3163 0.4492 0.7041 0.5186
22 4 0 1 A LFX -0.0052 0.0027 -1.9131 0.9443
23 4 0 1 A PRES -0.0654 0.0205 -3.1986 0.99865. Piecewise Constant Exponential Models
Book chapter 5, pp. 119-130. Command files ehe1 to
ehe3.
The rate is constant within periods and free between them, which buys duration dependence without committing to a shape.
Box 17The piecewise constant helper
pc <- function(f, data, model = "exponential_periods", tp = "0 (12) 96", ...)
tda_rate(f, data, model = model, options = list(tp = tp), ...)ehe1 is the null model, ehe2 the same with covariates.
Box 18ehe1, ehe2: piecewise constant exponential [book pp. 122-127, Boxes 5.2.1-5.2.2, 5.3.1-5.3.2]
ehe1 <- pc(Surv(TFP, DES) ~ 1, d, prate = "0(1)300")
ehe1Call: tda_rate(formula = f, data = data, model = model, prate = "0(1)300",
options = list(tp = tp))
Episodes: 600
Model: exponential_periods (TDA 3)
logLik (starting values): -2514.02
logLik: -2456.609
Converged in 5 iterations
Idx SN Org Des MT Variable Coeff Error C/Error Signif
1 1 0 1 A Period-1 -4.6829 0.1260 -37.1691 1.0000
2 1 0 1 A Period-2 -4.0420 0.1021 -39.6029 1.0000
3 1 0 1 A Period-3 -3.7854 0.1026 -36.8954 1.0000
4 1 0 1 A Period-4 -4.3178 0.1508 -28.6410 1.0000
5 1 0 1 A Period-5 -4.2955 0.1622 -26.4792 1.0000
6 1 0 1 A Period-6 -4.3690 0.1857 -23.5278 1.0000
7 1 0 1 A Period-7 -4.7274 0.2425 -19.4915 1.0000
8 1 0 1 A Period-8 -5.6864 0.4082 -13.9289 1.0000
9 1 0 1 A Period-9 -5.1650 0.1195 -43.2137 1.0000ehe1$episodes # the episodes TDA built sn org des episodes weighted mean.duration ts.min tf.max excluded
1 1 0 0 142 142 128.17606 0 428 FALSE
2 1 0 1 458 458 49.30349 0 350 FALSEehe1$periods # the period boundaries and the events in each start end starting ending events
1 0 12 600 76 63
2 12 24 0 104 96
3 24 36 0 108 95
4 36 48 0 55 44
5 48 60 0 44 38
6 60 72 0 38 29
7 72 84 0 22 17
8 84 96 0 9 6
9 96 Inf 0 144 70ehe2 <- pc(with_x("Surv(TFP, DES)"), d)
ehe2Call: tda_rate(formula = f, data = data, model = model, options = list(tp = tp))
Episodes: 600
Model: exponential_periods (TDA 3)
logLik (starting values): -2514.02
logLik: -2417.175
Converged in 6 iterations
Idx SN Org Des MT Variable Coeff Error C/Error Signif
1 1 0 1 A Period-1 -4.6053 0.3019 -15.2522 1.0000
2 1 0 1 A Period-2 -3.9236 0.2922 -13.4281 1.0000
3 1 0 1 A Period-3 -3.6302 0.2904 -12.5000 1.0000
4 1 0 1 A Period-4 -4.1514 0.3101 -13.3852 1.0000
5 1 0 1 A Period-5 -4.1106 0.3157 -13.0201 1.0000
6 1 0 1 A Period-6 -4.1584 0.3292 -12.6327 1.0000
7 1 0 1 A Period-7 -4.5272 0.3609 -12.5459 1.0000
8 1 0 1 A Period-8 -5.4741 0.4876 -11.2257 1.0000
9 1 0 1 A Period-9 -4.9094 0.2887 -17.0067 1.0000
10 1 0 1 A EDU 0.0683 0.0249 2.7409 0.9939
11 1 0 1 A COHO2 0.4841 0.1155 4.1931 1.0000
12 1 0 1 A COHO3 0.3823 0.1228 3.1146 0.9982
13 1 0 1 A LFX -0.0037 0.0009 -3.9943 0.9999
14 1 0 1 A PNOJ 0.0618 0.0442 1.3987 0.8381
15 1 0 1 A PRES -0.0266 0.0055 -4.8228 1.0000ehe2$periods start end starting ending events
1 0 12 600 76 63
2 12 24 0 104 96
3 24 36 0 108 95
4 36 48 0 55 44
5 48 60 0 44 38
6 60 72 0 38 29
7 72 84 0 22 17
8 84 96 0 9 6
9 96 Inf 0 144 70ehe1a plots the rate. The steps are the periods.
(ehe2p is not a second rate plot but the pseudoresidual
check of ehe2, which comes below with the rest of section
8.2.)
Box 19ehe1a: the estimated rate [book pp. 122-124, Boxes 5.2.1-5.2.2; Fig. 5.2.1]
er <- tda_rates(ehe1)
p <- tda_ps(er, width = 80, height = 40, xlim = c(0, 300), ylim = c(0, 0.03))
p <- tda_pl(p, "plxa", sc = 60, ic = 5)
p <- tda_pl(p, "plya", sc = 0.01, ic = 2)
p <- tda_pl_lines(p, "time", "rate")
plot(tda_pl_text(p, "Piecewise Constant Exponential Rate", at = c(60, 0.025)))
ehe3 uses three periods chosen by hand rather than a regular grid, which is model 16 rather than model 3.
Box 20ehe3: three periods [book pp. 129-130, Boxes 5.4.1-5.4.2]
ehe3 <- pc(with_x("Surv(TFP, DES)"), d,
model = "exponential_periods2", tp = "0,24,60")
ehe3Call: tda_rate(formula = f, data = data, model = model, options = list(tp = tp))
Episodes: 600
Model: exponential_periods2 (TDA 16)
logLik (starting values): -2514.02
logLik: -2434.896
Converged in 6 iterations
Idx SN Org Des MT P Variable Coeff Error C/Error Signif
1 1 0 1 A 1.0000 Period-1 -4.4627 0.4522 -9.8688 1.0000
2 1 0 1 A 2.0000 Period-2 -3.7235 0.4451 -8.3647 1.0000
3 1 0 1 A 3.0000 Period-3 -5.0902 0.5780 -8.8058 1.0000
4 1 0 1 A 1.0000 EDU 0.1016 0.0384 2.6455 0.9918
5 1 0 1 A 2.0000 EDU 0.0548 0.0435 1.2590 0.7920
6 1 0 1 A 3.0000 EDU 0.0470 0.0494 0.9508 0.6583
7 1 0 1 A 1.0000 COHO2 0.5849 0.2038 2.8692 0.9959
8 1 0 1 A 2.0000 COHO2 0.5112 0.1864 2.7426 0.9939
9 1 0 1 A 3.0000 COHO2 0.4116 0.2192 1.8776 0.9396
10 1 0 1 A 1.0000 COHO3 0.4666 0.2061 2.2638 0.9764
11 1 0 1 A 2.0000 COHO3 0.2824 0.1939 1.4562 0.8547
12 1 0 1 A 3.0000 COHO3 0.6422 0.2518 2.5505 0.9892
13 1 0 1 A 1.0000 LFX -0.0059 0.0017 -3.4044 0.9993
14 1 0 1 A 2.0000 LFX -0.0031 0.0014 -2.1951 0.9718
15 1 0 1 A 3.0000 LFX -0.0017 0.0018 -0.9437 0.6547
16 1 0 1 A 1.0000 PNOJ 0.1668 0.0742 2.2486 0.9755
17 1 0 1 A 2.0000 PNOJ 0.0326 0.0697 0.4683 0.3604
18 1 0 1 A 3.0000 PNOJ -0.0445 0.0917 -0.4856 0.3727
19 1 0 1 A 1.0000 PRES -0.0336 0.0089 -3.7938 0.9999
20 1 0 1 A 2.0000 PRES -0.0265 0.0091 -2.9246 0.9966
21 1 0 1 A 3.0000 PRES -0.0121 0.0113 -1.0723 0.7164# the two tables TDA prints before the estimates: the episodes it built,
# and how they fall across the periods
ehe3$episodes sn org des episodes weighted mean.duration ts.min tf.max excluded
1 1 0 0 142 142 128.17606 0 428 FALSE
2 1 0 1 458 458 49.30349 0 350 FALSEehe3$periods start end starting ending events
1 0 24 600 180 159
2 24 60 0 207 177
3 60 Inf 0 213 1226. Exponential Models with Time-dependent Covariates
Book chapter 6, pp. 131-157. Episode splitting; command files
ehf1 to ehf6.
ehf1 splits every episode at the date of marriage
and writes the pieces out (rrdat.s1, Box 6.4.2);
ehf5 splits on a fixed grid of 60 months
(rrdat.d60, Box 6.5.2). tda_split() is both,
at= for a per-episode split time and grid= for
the fixed one.
Box 21ehf1, ehf5: the split episode data [book pp. 143-154, Boxes 6.4.1-6.4.2, 6.5.1-6.5.2]
ds <- d
ds$MarrDate <- ifelse(ds$TMAR <= 0, 10000, ds$TMAR - ds$TStart)
sp <- tda_split(stats::reformulate(c(XVARS, "SEX"),
response = "Surv(TFP, DES)"),
ds, at = "MarrDate")
cat(sprintf("%d episodes split into %d\n", nrow(ds), nrow(sp)))600 episodes split into 761attr(sp, "episodes") # the episodes epdat reports before writing them sn org des episodes weighted mean.duration ts.min tf.max excluded
1 1 0 0 142 142 128.17606 0 428 FALSE
2 1 0 1 458 458 49.30349 0 350 FALSEhead(sp, 6) episode case subsample transition org des ts tf EDU COHO2 COHO3 LFX PNOJ PRES SEX MarrDate
1 1 1 2 1 0 0 0 124 17 0 0 0 0 34 1 124
2 1 1 2 2 0 0 124 428 17 0 0 0 0 34 1 124
3 2 1 1 1 0 1 0 46 10 0 0 0 0 22 2 169
4 3 1 1 1 0 1 0 34 10 0 0 46 1 46 2 123
5 4 1 2 1 0 0 0 89 10 0 0 80 2 46 2 89
6 4 1 2 2 0 1 89 220 10 0 0 80 2 46 2 89sp6 <- tda_split(with_x("Surv(TFP, DES)"), d, grid = seq(60, 420, by = 60))
cat(sprintf("%d episodes split into %d on the 60-month grid\n",
nrow(d), nrow(sp6)))600 episodes split into 1021 on the 60-month gridattr(sp6, "episodes") sn org des episodes weighted mean.duration ts.min tf.max excluded
1 1 0 0 142 142 128.17606 0 428 FALSE
2 1 0 1 458 458 49.30349 0 350 FALSEhead(sp6, 6) episode case subsample transition org des ts tf EDU COHO2 COHO3 LFX PNOJ PRES
1 1 1 8 1 0 0 0 60 17 0 0 0 0 34
2 1 1 8 2 0 0 60 120 17 0 0 0 0 34
3 1 1 8 3 0 0 120 180 17 0 0 0 0 34
4 1 1 8 4 0 0 180 240 17 0 0 0 0 34
5 1 1 8 5 0 0 240 300 17 0 0 0 0 34
6 1 1 8 6 0 0 300 360 17 0 0 0 0 34A covariate that changes during an episode is handled by cutting the
episode where it changes. Marriage is the example: MarrDate
is the month it happens, or past the end of the episode if it never
does. ehf2, ehf3,
ehf4 and ehf6 are the models on the
pieces.
Box 22ehf2, ehf3: episode splitting at marriage [book pp. 145-149, Boxes 6.4.3-6.4.6]
sp$MARR <- as.integer(sp$MarrDate <= sp$ts)
cat(sprintf("%d episodes split into %d pieces\n", nrow(ds), nrow(sp)))600 episodes split into 761 pieces# Surv(ts, tf, des), not Surv(tf - ts, des): a piece that begins at ts is
# left truncated, and a model that is told otherwise puts it into risk
# sets it was never in. An exponential rate only ever sees the difference
# and so is unharmed, but the Cox models of ehi3 to ehi5 below are not.
ehf2 <- tda_rate(Surv(ts, tf, des) ~ MARR, sp)
ehf2$episodes sn org des episodes weighted mean.duration ts.min tf.max excluded
1 1 0 0 303 303 76.35313 0 428 FALSE
2 1 0 1 458 458 38.53057 0 350 FALSEehf2Call: tda_rate(formula = Surv(ts, tf, des) ~ MARR, data = sp)
Episodes: 761
Model: exponential (TDA 2)
logLik (starting values): -2514.02
logLik: -2498.534
Converged in 4 iterations
Idx SN Org Des MT Variable Coeff Error C/Error Signif
1 1 0 1 A Constant -4.2129 0.0639 -65.9430 1.0000
2 1 0 1 A MARR -0.5212 0.0937 -5.5629 1.0000ehf3 <- tda_rate(stats::reformulate(c("MARR", XVARS),
response = "Surv(ts, tf, des)"), sp)
ehf3$episodes sn org des episodes weighted mean.duration ts.min tf.max excluded
1 1 0 0 303 303 76.35313 0 428 FALSE
2 1 0 1 458 458 38.53057 0 350 FALSEehf3Call: tda_rate(formula = stats::reformulate(c("MARR", XVARS), response = "Surv(ts, tf, des)"),
data = sp)
Episodes: 761
Model: exponential (TDA 2)
logLik (starting values): -2514.02
logLik: -2460.248
Converged in 5 iterations
Idx SN Org Des MT Variable Coeff Error C/Error Signif
1 1 0 1 A Constant -4.3878 0.2756 -15.9222 1.0000
2 1 0 1 A MARR -0.3447 0.1022 -3.3716 0.9993
3 1 0 1 A EDU 0.0770 0.0246 3.1358 0.9983
4 1 0 1 A COHO2 0.5971 0.1137 5.2534 1.0000
5 1 0 1 A COHO3 0.6138 0.1188 5.1651 1.0000
6 1 0 1 A LFX -0.0023 0.0010 -2.3731 0.9824
7 1 0 1 A PNOJ 0.0645 0.0440 1.4674 0.8577
8 1 0 1 A PRES -0.0274 0.0055 -4.9533 1.0000ehf4 distinguishes married men from married women.
Box 23ehf4: marriage, by sex [book pp. 150-151, Boxes 6.4.7-6.4.8]
sp$MarrMen <- as.integer(sp$SEX == 1 & sp$MarrDate <= sp$ts)
ehf4 <- tda_rate(stats::reformulate(c(XVARS, "MARR", "MarrMen"),
response = "Surv(ts, tf, des)"), sp)
ehf4Call: tda_rate(formula = stats::reformulate(c(XVARS, "MARR", "MarrMen"),
response = "Surv(ts, tf, des)"), data = sp)
Episodes: 761
Model: exponential (TDA 2)
logLik (starting values): -2514.02
logLik: -2434.175
Converged in 6 iterations
Idx SN Org Des MT Variable Coeff Error C/Error Signif
1 1 0 1 A Constant -4.5842 0.2694 -17.0163 1.0000
2 1 0 1 A EDU 0.0943 0.0234 4.0227 0.9999
3 1 0 1 A COHO2 0.5423 0.1132 4.7916 1.0000
4 1 0 1 A COHO3 0.5304 0.1192 4.4494 1.0000
5 1 0 1 A LFX -0.0031 0.0010 -3.1550 0.9984
6 1 0 1 A PNOJ 0.1029 0.0445 2.3126 0.9793
7 1 0 1 A PRES -0.0264 0.0054 -4.9363 1.0000
8 1 0 1 A MARR 0.2735 0.1231 2.2209 0.9736
9 1 0 1 A MarrMen -1.0231 0.1401 -7.3019 1.0000ehf4$episodes sn org des episodes weighted mean.duration ts.min tf.max excluded
1 1 0 0 303 303 76.35313 0 428 FALSE
2 1 0 1 458 458 38.53057 0 350 FALSEehf6 is the model on the episodes
ehf5 split on a fixed grid – at 60, 120, … 420 months –
so that a covariate measured in calendar time can change within a job.
grid= takes those times; at= is for a split
time that differs per episode, as marriage does. LFX60 is
the sub-episode’s starting time, which is labour force experience
accumulated within the job.
Box 24ehf6: splitting on a fixed grid [book pp. 155-156, Boxes 6.5.3-6.5.4]
sp6$LFX60 <- sp6$ts
cat(sprintf("%d episodes split into %d pieces\n", nrow(d), nrow(sp6)))600 episodes split into 1021 piecesehf6 <- tda_rate(stats::reformulate(c(XVARS, "LFX60"),
response = "Surv(ts, tf, des)"), sp6)
ehf6Call: tda_rate(formula = stats::reformulate(c(XVARS, "LFX60"), response = "Surv(ts, tf, des)"),
data = sp6)
Episodes: 1021
Model: exponential (TDA 2)
logLik (starting values): -2514.02
logLik: -2432.308
Converged in 6 iterations
Idx SN Org Des MT Variable Coeff Error C/Error Signif
1 1 0 1 A Constant -4.0138 0.2764 -14.5194 1.0000
2 1 0 1 A EDU 0.0646 0.0248 2.6051 0.9908
3 1 0 1 A COHO2 0.4073 0.1149 3.5462 0.9996
4 1 0 1 A COHO3 0.2946 0.1216 2.4224 0.9846
5 1 0 1 A LFX -0.0039 0.0009 -4.2429 1.0000
6 1 0 1 A PNOJ 0.0647 0.0440 1.4705 0.8586
7 1 0 1 A PRES -0.0253 0.0054 -4.6474 1.0000
8 1 0 1 A LFX60 -0.0071 0.0010 -7.1098 1.0000ehf6$episodes # the split episodes, as TDA reports them sn org des episodes weighted mean.duration ts.min tf.max excluded
1 1 0 0 563 563 51.61812 0 428 FALSE
2 1 0 1 458 458 25.59170 0 350 FALSEBox 6.5.5 of the book does the same with finer grids – 24 and 12 months – and shows how little the estimates move with the grid.
Box 25Box 6.5.5: the grid at 60, 24 and 12 months [book p. 157, Box 6.5.5]
grid_fit <- function(step) {
spg <- tda_split(with_x("Surv(TFP, DES)"), d, grid = seq(step, 420, by = step))
spg$LFXG <- spg$ts
fit <- tda_rate(stats::reformulate(c(XVARS, "LFXG"),
response = "Surv(ts, tf, des)"), spg)
e <- fit$estimates
stats::setNames(data.frame(e$Coeff, e$`C/Error`), c("Coeff", "T-Stat"))
}
# the three grids side by side, coefficient and t statistic for each
g <- lapply(c(60, 24, 12), grid_fit)
out <- do.call(cbind, g)
names(out) <- paste(rep(c("60", "24", "12"), each = 2), names(out))
rownames(out) <- c("Const", XVARS, "LFXG")
print(format(out, digits = 5)) 60 Coeff 60 T-Stat 24 Coeff 24 T-Stat 12 Coeff 12 T-Stat
Const -4.0137800 -14.5194 -3.9498113 -14.1626 -3.9187212 -14.0056
EDU 0.0646116 2.6051 0.0628974 2.5367 0.0626471 2.5264
COHO2 0.4073219 3.5462 0.4059659 3.5276 0.4049094 3.5177
COHO3 0.2946223 2.4224 0.3006367 2.4641 0.2996529 2.4555
LFX -0.0039428 -4.2429 -0.0038975 -4.1976 -0.0038894 -4.1903
PNOJ 0.0647431 1.4705 0.0632658 1.4359 0.0628425 1.4263
PRES -0.0253234 -4.6474 -0.0249290 -4.5755 -0.0248351 -4.5595
LFXG -0.0070767 -7.1098 -0.0063135 -6.9933 -0.0062674 -7.01407. Parametric Models of Time-dependence
Book chapter 7, pp. 159-212. Gompertz, Weibull, log-logistic,
log-normal and sickle; command files ehg1 to
ehg13.
Five families, each fitted without covariates and then with them.
on= selects which parameter of the distribution the
covariates attach to — xa, xb or
xc — which is the difference between the
ehg2/ehg3/ ehg4 group of
files.
Box 26The two fits each family gets
fit1 <- function(model, ...)
tda_rate(Surv(TFP, DES) ~ 1, d, model = model, prate = "0(1)300", ...)
fitx <- function(model, on = "xa", ...)
tda_rate(with_x("Surv(TFP, DES)"), d, model = model, on = on, ...)ehg1, ehg5, ehg8, ehg10 and ehg12 are the null models: Gompertz, Weibull, log-logistic, log-normal and sickle.
Box 27ehg1, ehg5, ehg8, ehg10, ehg12: five families [book pp. 183-210, Boxes 7.2.1-7.2.2, 7.3.1-7.3.2, 7.4.1-7.4.2, 7.5.1-7.5.2, 7.6.1-7.6.2]
models <- c("gompertz", "weibull", "loglogistic", "lognormal", "sickle")
fits <- lapply(models, fit1)
names(fits) <- models
for (m in models) print(fits[[m]])Call: tda_rate(formula = Surv(TFP, DES) ~ 1, data = d, model = model,
prate = "0(1)300")
Episodes: 600
Model: gompertz (TDA 6)
logLik (starting values): -2557.365
logLik: -2474.506
Converged in 5 iterations
Idx SN Org Des MT Variable Coeff Error C/Error Signif
1 1 0 1 B Constant -4.0729 0.0633 -64.3172 1.0000
2 1 0 1 C Constant -0.0067 0.0009 -7.7201 1.0000
Call: tda_rate(formula = Surv(TFP, DES) ~ 1, data = d, model = model,
prate = "0(1)300")
Episodes: 600
Model: weibull (TDA 7)
logLik (starting values): -2514.02
logLik: -2504.899
Converged in 4 iterations
Idx SN Org Des MT Variable Coeff Error C/Error Signif
1 1 0 1 A Constant -4.4616 0.0544 -82.0244 1.0000
2 1 0 1 B Constant -0.1477 0.0360 -4.1025 1.0000
Call: tda_rate(formula = Surv(TFP, DES) ~ 1, data = d, model = model,
prate = "0(1)300")
Episodes: 600
Model: loglogistic (TDA 9)
logLik (starting values): -2521.785
logLik: -2460.487
Converged in 5 iterations
Idx SN Org Des MT Variable Coeff Error C/Error Signif
1 1 0 1 A Constant -3.8434 0.0549 -70.0613 1.0000
2 1 0 1 B Constant 0.2918 0.0385 7.5742 1.0000
Call: tda_rate(formula = Surv(TFP, DES) ~ 1, data = d, model = model,
prate = "0(1)300")
Episodes: 600
Model: lognormal (TDA 12)
logLik (starting values): -2603.168
logLik: -2456.138
Converged in 6 iterations
Idx SN Org Des MT Variable Coeff Error C/Error Signif
1 1 0 1 A Constant 3.8852 0.0548 70.8561 1.0000
2 1 0 1 B Constant 0.2436 0.0342 7.1252 1.0000
Call: tda_rate(formula = Surv(TFP, DES) ~ 1, data = d, model = model,
prate = "0(1)300")
Episodes: 600
Model: sickle (TDA 8)
logLik (starting values): -2897.762
logLik: -2486.366
Converged in 6 iterations
Idx SN Org Des MT Variable Coeff Error C/Error Signif
1 1 0 1 A Constant -6.5843 0.0775 -84.9960 1.0000
2 1 0 1 B Constant 3.6467 0.0481 75.8802 1.0000# all five are fitted to the same episodes, which TDA reports before
# each set of estimates
fits$gompertz$episodes sn org des episodes weighted mean.duration ts.min tf.max excluded
1 1 0 0 142 142 128.17606 0 428 FALSE
2 1 0 1 458 458 49.30349 0 350 FALSEehg1a, ehg5a and the rest plot each fitted rate, which is where the families differ from each other and from the exponential. All five command files set the same axes – 0 to 300 months against a rate of 0 to 0.02 – so the shapes can be compared directly. The book prints them one to a section, beside the model they belong to: Figures 7.2.2, 7.3.2, 7.4.2, 7.5.2 and 7.6.2, in that order. The exponential’s constant 0.0112 would be a horizontal line across the middle of each.
Box 28ehg1a, ehg5a, ehg8a, ehg10a, ehg12a: the five fitted rates [book pp. 185-211, Fig. 7.2.2, Fig. 7.3.2, Fig. 7.4.2, Fig. 7.5.2, Fig. 7.6.2]
p <- tda_ps(tda_rates(fits$gompertz), file = "ehg1a.ps", width = 80, height = 40,
xlim = c(0, 300), ylim = c(0, 0.02))
p <- tda_pl(p, "plxa", sc = 60, ic = 5)
p <- tda_pl(p, "plya", sc = 0.01, ic = 5)
plot(tda_pl_lines(p, "time", "rate"))
p <- tda_ps(tda_rates(fits$weibull), file = "ehg5a.ps", width = 80, height = 40,
xlim = c(0, 300), ylim = c(0, 0.02))
p <- tda_pl(p, "plxa", sc = 60, ic = 5)
p <- tda_pl(p, "plya", sc = 0.01, ic = 5)
plot(tda_pl_lines(p, "time", "rate"))
p <- tda_ps(tda_rates(fits$loglogistic), file = "ehg8a.ps", width = 80, height = 40,
xlim = c(0, 300), ylim = c(0, 0.02))
p <- tda_pl(p, "plxa", sc = 60, ic = 5)
p <- tda_pl(p, "plya", sc = 0.01, ic = 5)
plot(tda_pl_lines(p, "time", "rate"))
p <- tda_ps(tda_rates(fits$lognormal), file = "ehg10a.ps", width = 80, height = 40,
xlim = c(0, 300), ylim = c(0, 0.02))
p <- tda_pl(p, "plxa", sc = 60, ic = 5)
p <- tda_pl(p, "plya", sc = 0.01, ic = 5)
plot(tda_pl_lines(p, "time", "rate"))
p <- tda_ps(tda_rates(fits$sickle), file = "ehg12a.ps", width = 80, height = 40,
xlim = c(0, 300), ylim = c(0, 0.02))
p <- tda_pl(p, "plxa", sc = 60, ic = 5)
p <- tda_pl(p, "plya", sc = 0.01, ic = 5)
plot(tda_pl_lines(p, "time", "rate"))
ehg2, ehg6, ehg9,
ehg11 and ehg13 add the covariates. In
the Gompertz-Makeham model of ehg2 they go on
xb, the term the rate’s growth depends on, not on the
Makeham constant xa; the other four families put them on
xa.
Box 29ehg2, ehg6, ehg9, ehg11, ehg13: with covariates [book pp. 186-212, Boxes 7.2.3-7.2.4, 7.3.3-7.3.4, 7.4.3-7.4.4, 7.5.3-7.5.4, 7.6.3-7.6.4]
sickle13 <- fitx("sickle")
sickle13$episodes # the same episodes for all five, as TDA reports sn org des episodes weighted mean.duration ts.min tf.max excluded
1 1 0 0 142 142 128.17606 0 428 FALSE
2 1 0 1 458 458 49.30349 0 350 FALSEfitx("gompertz", on = "xb", prate = "0(10)300")Call: tda_rate(formula = with_x("Surv(TFP, DES)"), data = d, model = model,
on = on, prate = "0(10)300")
Episodes: 600
Model: gompertz (TDA 6)
logLik (starting values): -2557.365
logLik: -2437.975
Converged in 5 iterations
Idx SN Org Des MT Variable Coeff Error C/Error Signif
1 1 0 1 B Constant -3.9105 0.2814 -13.8985 1.0000
2 1 0 1 B EDU 0.0634 0.0248 2.5565 0.9894
3 1 0 1 B COHO2 0.4130 0.1152 3.5849 0.9997
4 1 0 1 B COHO3 0.3119 0.1222 2.5512 0.9893
5 1 0 1 B LFX -0.0039 0.0009 -4.1613 1.0000
6 1 0 1 B PNOJ 0.0627 0.0441 1.4241 0.8456
7 1 0 1 B PRES -0.0250 0.0055 -4.5893 1.0000
8 1 0 1 C Constant -0.0059 0.0009 -6.7272 1.0000fitx("weibull", prate = "0(10)300")Call: tda_rate(formula = with_x("Surv(TFP, DES)"), data = d, model = model,
on = on, prate = "0(10)300")
Episodes: 600
Model: weibull (TDA 7)
logLik (starting values): -2514.02
logLik: -2462.759
Converged in 6 iterations
Idx SN Org Des MT Variable Coeff Error C/Error Signif
1 1 0 1 A Constant -4.4063 0.3060 -14.3996 1.0000
2 1 0 1 A EDU 0.0779 0.0270 2.8848 0.9961
3 1 0 1 A COHO2 0.6062 0.1242 4.8821 1.0000
4 1 0 1 A COHO3 0.5777 0.1303 4.4318 1.0000
5 1 0 1 A LFX -0.0036 0.0010 -3.5040 0.9995
6 1 0 1 A PNOJ 0.0637 0.0484 1.3184 0.8126
7 1 0 1 A PRES -0.0292 0.0060 -4.8396 1.0000
8 1 0 1 B Constant -0.0903 0.0364 -2.4830 0.9870fitx("loglogistic")Call: tda_rate(formula = with_x("Surv(TFP, DES)"), data = d, model = model,
on = on)
Episodes: 600
Model: loglogistic (TDA 9)
logLik (starting values): -2521.785
logLik: -2418.904
Converged in 5 iterations
Idx SN Org Des MT Variable Coeff Error C/Error Signif
1 1 0 1 A Constant -3.7785 0.2926 -12.9144 1.0000
2 1 0 1 A EDU 0.0819 0.0269 3.0485 0.9977
3 1 0 1 A COHO2 0.5338 0.1249 4.2726 1.0000
4 1 0 1 A COHO3 0.4005 0.1318 3.0382 0.9976
5 1 0 1 A LFX -0.0043 0.0009 -4.6566 1.0000
6 1 0 1 A PNOJ 0.0970 0.0476 2.0368 0.9583
7 1 0 1 A PRES -0.0297 0.0056 -5.2900 1.0000
8 1 0 1 B Constant 0.3615 0.0386 9.3736 1.0000fitx("lognormal")Call: tda_rate(formula = with_x("Surv(TFP, DES)"), data = d, model = model,
on = on)
Episodes: 600
Model: lognormal (TDA 12)
logLik (starting values): -2603.168
logLik: -2415.882
Converged in 6 iterations
Idx SN Org Des MT Variable Coeff Error C/Error Signif
1 1 0 1 A Constant 3.8567 0.2916 13.2280 1.0000
2 1 0 1 A EDU -0.0813 0.0270 -3.0127 0.9974
3 1 0 1 A COHO2 -0.5164 0.1247 -4.1422 1.0000
4 1 0 1 A COHO3 -0.4514 0.1295 -3.4868 0.9995
5 1 0 1 A LFX 0.0039 0.0009 4.4413 1.0000
6 1 0 1 A PNOJ -0.0849 0.0464 -1.8296 0.9327
7 1 0 1 A PRES 0.0287 0.0056 5.1330 1.0000
8 1 0 1 B Constant 0.1794 0.0342 5.2437 1.0000fitx("sickle")Call: tda_rate(formula = with_x("Surv(TFP, DES)"), data = d, model = model,
on = on)
Episodes: 600
Model: sickle (TDA 8)
logLik (starting values): -2897.762
logLik: -2443.935
Converged in 8 iterations
Idx SN Org Des MT Variable Coeff Error C/Error Signif
1 1 0 1 A Constant -6.4245 0.2831 -22.6909 1.0000
2 1 0 1 A EDU 0.0767 0.0250 3.0711 0.9979
3 1 0 1 A COHO2 0.4259 0.1149 3.7073 0.9998
4 1 0 1 A COHO3 0.2785 0.1206 2.3095 0.9791
5 1 0 1 A LFX -0.0043 0.0009 -4.5725 1.0000
6 1 0 1 A PNOJ 0.0782 0.0438 1.7867 0.9260
7 1 0 1 A PRES -0.0288 0.0055 -5.1914 1.0000
8 1 0 1 B Constant 3.6897 0.0502 73.4793 1.0000ehg3 and ehg4 put the same
covariates on the other parameters as well: xb and
xc, then all three. ehg7 does the same for
the Weibull, on xa and xb.
Box 30ehg3, ehg4, ehg7: covariates on more than one term [book pp. 188-197, Boxes 7.2.5-7.2.8, 7.3.5-7.3.6]
ehg4 <- fitx("gompertz", on = c("xa", "xb", "xc"))
ehg4$episodes # the same episodes for every model in this group sn org des episodes weighted mean.duration ts.min tf.max excluded
1 1 0 0 142 142 128.17606 0 428 FALSE
2 1 0 1 458 458 49.30349 0 350 FALSEfitx("gompertz", on = c("xb", "xc"))Call: tda_rate(formula = with_x("Surv(TFP, DES)"), data = d, model = model,
on = on)
Episodes: 600
Model: gompertz (TDA 6)
logLik (starting values): -2557.365
logLik: -2433.387
Converged in 6 iterations
Idx SN Org Des MT Variable Coeff Error C/Error Signif
1 1 0 1 B Constant -4.0167 0.3663 -10.9652 1.0000
2 1 0 1 B EDU 0.0918 0.0318 2.8843 0.9961
3 1 0 1 B COHO2 0.5074 0.1547 3.2789 0.9990
4 1 0 1 B COHO3 0.2069 0.1687 1.2269 0.7801
5 1 0 1 B LFX -0.0051 0.0013 -3.9165 0.9999
6 1 0 1 B PNOJ 0.1513 0.0604 2.5060 0.9878
7 1 0 1 B PRES -0.0323 0.0071 -4.5315 1.0000
8 1 0 1 C Constant -0.0054 0.0057 -0.9497 0.6577
9 1 0 1 C EDU -0.0005 0.0005 -1.1416 0.7464
10 1 0 1 C COHO2 -0.0019 0.0022 -0.8685 0.6149
11 1 0 1 C COHO3 0.0030 0.0028 1.0626 0.7120
12 1 0 1 C LFX 0.0000 0.0000 1.3730 0.8303
13 1 0 1 C PNOJ -0.0019 0.0010 -2.0254 0.9572
14 1 0 1 C PRES 0.0002 0.0001 1.6056 0.8916fitx("gompertz", on = c("xa", "xb", "xc"))Call: tda_rate(formula = with_x("Surv(TFP, DES)"), data = d, model = model,
on = on)
Episodes: 600
Model: gompertz (TDA 6)
logLik (starting values): -2538.685
logLik: -2423.848
Converged in 18 iterations
Idx SN Org Des MT Variable Coeff Error C/Error Signif
1 1 0 1 A Constant -14.9029 3.6611 -4.0706 1.0000
2 1 0 1 A EDU -0.1348 0.2272 -0.5932 0.4470
3 1 0 1 A COHO2 2.2350 1.7415 1.2834 0.8006
4 1 0 1 A COHO3 -3.4706 14.2396 -0.2437 0.1926
5 1 0 1 A LFX 0.0162 0.0082 1.9770 0.9520
6 1 0 1 A PNOJ -2.9917 0.8958 -3.3397 0.9992
7 1 0 1 A PRES 0.2087 0.0601 3.4714 0.9995
8 1 0 1 B Constant -3.8740 0.3855 -10.0491 1.0000
9 1 0 1 B EDU 0.0988 0.0334 2.9581 0.9969
10 1 0 1 B COHO2 0.4783 0.1673 2.8581 0.9957
11 1 0 1 B COHO3 0.2640 0.1741 1.5164 0.8706
12 1 0 1 B LFX -0.0057 0.0014 -3.9031 0.9999
13 1 0 1 B PNOJ 0.1898 0.0624 3.0432 0.9977
14 1 0 1 B PRES -0.0407 0.0080 -5.1086 1.0000
15 1 0 1 C Constant -0.0045 0.0064 -0.7146 0.5251
16 1 0 1 C EDU -0.0005 0.0005 -1.0523 0.7073
17 1 0 1 C COHO2 -0.0024 0.0024 -0.9801 0.6729
18 1 0 1 C COHO3 0.0031 0.0029 1.0850 0.7221
19 1 0 1 C LFX 0.0000 0.0000 1.4054 0.8401
20 1 0 1 C PNOJ -0.0018 0.0010 -1.8275 0.9324
21 1 0 1 C PRES 0.0002 0.0001 1.2268 0.7801fitx("weibull", on = c("xa", "xb"))Call: tda_rate(formula = with_x("Surv(TFP, DES)"), data = d, model = model,
on = on)
Episodes: 600
Model: weibull (TDA 7)
logLik (starting values): -2514.02
logLik: -2459.595
Converged in 7 iterations
Idx SN Org Des MT Variable Coeff Error C/Error Signif
1 1 0 1 A Constant -4.4951 0.3073 -14.6293 1.0000
2 1 0 1 A EDU 0.0858 0.0274 3.1351 0.9983
3 1 0 1 A COHO2 0.5985 0.1285 4.6586 1.0000
4 1 0 1 A COHO3 0.5582 0.1297 4.3040 1.0000
5 1 0 1 A LFX -0.0032 0.0011 -2.9481 0.9968
6 1 0 1 A PNOJ 0.0571 0.0496 1.1512 0.7503
7 1 0 1 A PRES -0.0290 0.0061 -4.7176 1.0000
8 1 0 1 B Constant -0.1613 0.2224 -0.7252 0.5317
9 1 0 1 B EDU -0.0027 0.0188 -0.1441 0.1146
10 1 0 1 B COHO2 0.0835 0.0856 0.9751 0.6705
11 1 0 1 B COHO3 0.2017 0.0914 2.2051 0.9726
12 1 0 1 B LFX 0.0010 0.0007 1.4443 0.8513
13 1 0 1 B PNOJ -0.0385 0.0363 -1.0614 0.7115
14 1 0 1 B PRES 0.0005 0.0041 0.1249 0.09948. Methods to Check Parametric Assumptions
Book chapter 8, pp. 213-228. The transformed survivor function
and the pseudoresiduals; command files ehh1 to
ehh4 and the *s/*p
pairs.
Each of these transforms the survivor function so that one distribution appears as a straight line, and regresses it on time. The closer to a line, the better that distribution fits.
The regression runs on the survivor function at five decimals, the
precision ehc5.ple carries: the command files read that
file rather than the estimate in memory, and the transformations are
sensitive enough for it to show in the t statistics.
The axes are each command file’s own.
Box 31ehh1 to ehh4: the four transformations [book pp. 215-218, Boxes 8.1.1-8.1.3]
s <- tda_survivor(pl, conf.int = NULL)
s <- s[s$survivor > 0 & s$survivor < 1 & s$time > 0, ]
s$survivor <- round(s$survivor, 5) # the precision ehc5.ple carries
checks <- list(
Exponential = list(y = log(s$survivor), x = s$time,
xlim = c(0, 300), ylim = c(-3, 0), xsc = 60, ysc = 1),
Weibull = list(y = log(-log(s$survivor)), x = log(s$time),
xlim = c(0, 6), ylim = c(-6, 2), xsc = 1, ysc = 1),
`Log-Logistic` = list(y = log(1 / s$survivor - 1), x = log(s$time),
xlim = c(0, 6), ylim = c(-6, 2), xsc = 1, ysc = 1),
`Log-Normal` = list(y = stats::qnorm(1 - s$survivor), x = log(s$time),
xlim = c(0, 6), ylim = c(-3, 2), xsc = 1, ysc = 1))
ehh <- lapply(checks, function(z) {
keep <- is.finite(z$y) & is.finite(z$x)
tda_lsreg(y ~ x, data.frame(x = z$x[keep], y = z$y[keep]))
})
# each regression in full: standard errors and t statistics with the
# two coefficients
for (nm in names(ehh)) {
cat("\n", nm, "\n", sep = "")
print(ehh[[nm]])
}
Exponential
Call: tda_lsreg(formula = y ~ x, data = data.frame(x = z$x[keep], y = z$y[keep]))
Cases: 130
Idx Wave Variable Coeff Error C/Error Signif
1 NA Intercept -0.3527 0.0280 -12.6000 1.0000
2 1.0000 x -0.0075 0.0002 -30.7721 1.0000
Weibull
Call: tda_lsreg(formula = y ~ x, data = data.frame(x = z$x[keep], y = z$y[keep]))
Cases: 130
Idx Wave Variable Coeff Error C/Error Signif
1 NA Intercept -4.4238 0.1135 -38.9803 1.0000
2 1.0000 x 1.0080 0.0270 37.3474 1.0000
Log-Logistic
Call: tda_lsreg(formula = y ~ x, data = data.frame(x = z$x[keep], y = z$y[keep]))
Cases: 130
Idx Wave Variable Coeff Error C/Error Signif
1 NA Intercept -5.1604 0.0771 -66.9338 1.0000
2 1.0000 x 1.3268 0.0183 72.3629 1.0000
Log-Normal
Call: tda_lsreg(formula = y ~ x, data = data.frame(x = z$x[keep], y = z$y[keep]))
Cases: 130
Idx Wave Variable Coeff Error C/Error Signif
1 NA Intercept -2.9220 0.0301 -97.2089 1.0000
2 1.0000 x 0.7578 0.0071 105.9948 1.0000The slope of the exponential line is the rate the transformation
implies, 0.0075, against the maximum likelihood estimate of
ehd1. Each command file ends by drawing its own fitted line
over the curve, the endpoints written into plotp by hand
(plotp = 0, -0.3527, 300, -2.6027 in ehh1.cf);
here the line comes from the fit above, and the four panels are laid out
as the book’s Figure 8.1.1.
Box 32ehh1, ehh2, ehh3, ehh4: the transformed survivor functions [book pp. 215-217, Boxes 8.1.1-8.1.2; Fig. 8.1.1]
z <- checks$`Exponential`
b <- coef(ehh$`Exponential`)
p <- tda_ps(data.frame(x = z$x, y = z$y), file = "ehh1.ps",
width = 80, height = 40, xlim = c(0, 300), ylim = c(-3, 0))
p <- tda_pl(p, "plxa", sc = 60, ic = 5)
p <- tda_pl(p, "plya", sc = 1, ic = 10)
p <- tda_pl(p, "plframe")
p <- tda_pl_lines(p, "x", "y")
p <- tda_pl_polyline(p, c(0, 300), b[[1L]] + b[[2L]] * c(0, 300))
ehh1 <- tda_pl_text(p, "Graphical Test (Exponential)", at = c(150.0, -0.30000000000000004))
z <- checks$`Weibull`
b <- coef(ehh$`Weibull`)
p <- tda_ps(data.frame(x = z$x, y = z$y), file = "ehh2.ps",
width = 80, height = 40, xlim = c(0, 6), ylim = c(-6, 2))
p <- tda_pl(p, "plxa", sc = 1, ic = 5)
p <- tda_pl(p, "plya", sc = 1, ic = 2)
p <- tda_pl(p, "plframe")
p <- tda_pl_lines(p, "x", "y")
p <- tda_pl_polyline(p, c(0, 6), b[[1L]] + b[[2L]] * c(0, 6))
ehh2 <- tda_pl_text(p, "Graphical Test (Weibull)", at = c(3.0, 1.2))
z <- checks$`Log-Logistic`
b <- coef(ehh$`Log-Logistic`)
p <- tda_ps(data.frame(x = z$x, y = z$y), file = "ehh3.ps",
width = 80, height = 40, xlim = c(0, 6), ylim = c(-6, 2))
p <- tda_pl(p, "plxa", sc = 1, ic = 5)
p <- tda_pl(p, "plya", sc = 1, ic = 2)
p <- tda_pl(p, "plframe")
p <- tda_pl_lines(p, "x", "y")
p <- tda_pl_polyline(p, c(0, 6), b[[1L]] + b[[2L]] * c(0, 6))
ehh3 <- tda_pl_text(p, "Graphical Test (Log-Logistic)", at = c(3.0, 1.2))
z <- checks$`Log-Normal`
b <- coef(ehh$`Log-Normal`)
p <- tda_ps(data.frame(x = z$x, y = z$y), file = "ehh4.ps",
width = 80, height = 40, xlim = c(0, 6), ylim = c(-3, 2))
p <- tda_pl(p, "plxa", sc = 1, ic = 5)
p <- tda_pl(p, "plya", sc = 1, ic = 2)
p <- tda_pl(p, "plframe")
p <- tda_pl_lines(p, "x", "y")
p <- tda_pl_polyline(p, c(0, 6), b[[1L]] + b[[2L]] * c(0, 6))
ehh4 <- tda_pl_text(p, "Graphical Test (Log-Normal)", at = c(3.0, 1.5))
r <- tda_combine_ps(list(list(ehh1, ehh2), list(ehh3, ehh4)), width = 170, height = 90)
tda_plot_ps(tda_read_ps(file.path(r$dir, "combined.ps")))
The book’s section 8.2 checks a model through its pseudoresiduals, in
three steps: residuals = TRUE writes them
(ehd2.res, Box 8.2.3); a product-limit estimate treats them
as durations (ehd2s.cf, Box 8.2.5); and
ehd2p.cf plots log S against the residual. If
the model is right the residuals are standard exponential, so that curve
should lie on a straight line of slope -1 through the origin.
ehd2p.cf draws that line beside it,
plotp = 0,0,3,-3; the plot is read as the distance between
the two.
Box 33ehd2.res, ehd2s: the pseudoresiduals and their survivor function [book pp. 97-225, Boxes 4.1.4-4.1.5, 8.2.1, 8.2.3-8.2.5]
ehd2r <- tda_rate(with_x("Surv(TFP, DES)"), d, residuals = TRUE)
head(ehd2r$residuals, 10) Case Org Des TS TF Rate Function Residual Weight
1 1 0 0 0 428 0.016122457 0.001007371 6.9004117 1
2 2 0 1 0 46 0.013133700 0.546538678 0.6041502 1
3 3 0 1 0 34 0.006149795 0.811319744 0.2090930 1
4 4 0 1 0 220 0.005858903 0.275557591 1.2889586 1
5 5 0 1 0 12 0.015308310 0.832185660 0.1836997 1
6 6 0 1 0 30 0.015640803 0.625487409 0.4692241 1
7 7 0 1 0 12 0.013875494 0.846617799 0.1665059 1
8 8 0 1 0 75 0.014176867 0.345326480 1.0632650 1
9 9 0 1 0 12 0.011855660 0.867388832 0.1422679 1
10 10 0 1 0 55 0.012105645 0.513856882 0.6658105 1tail(ehd2r$residuals, 6) Case Org Des TS TF Rate Function Residual Weight
595 595 0 1 0 119 0.008084062 0.3821266 0.9620033 1
596 596 0 1 0 66 0.005404227 0.6999972 0.3566790 1
597 597 0 0 0 112 0.005058187 0.5674986 0.5665169 1
598 598 0 0 0 103 0.017055411 0.1726123 1.7567073 1
599 599 0 1 0 22 0.019696756 0.6483474 0.4333286 1
600 600 0 1 0 13 0.019126466 0.7798575 0.2486441 1rp <- tda_ple(Surv(Residual, Des) ~ 1, ehd2r$residuals)
rp$episodes # the residuals treated as durations sn org des episodes weighted mean.duration ts.min tf.max excluded
1 1 0 0 142 142 1.2893610 0 6.900412 FALSE
2 1 0 1 458 458 0.6002418 0 3.595082 FALSEprint(rp, n = 10)Call: tda_ple(formula = Surv(Residual, Des) ~ 1, data = ehd2r$residuals)
SN 1. Transition: 0,1 - Product-Limit Estimation
id index time events censored n.risk survivor std.err cum.rate
0 0 0.00000000 0 0 600 1.0000000 0.000000000 0.000000000
0 1 0.01860143 1 1 599 0.9983306 0.001668055 0.001670844
0 2 0.01961462 1 0 598 0.9966611 0.002357013 0.003344485
0 3 0.03639071 1 2 595 0.9949860 0.002887553 0.005026571
0 4 0.03860349 1 0 594 0.9933110 0.003333319 0.006711491
0 5 0.04352060 1 1 592 0.9916331 0.003726134 0.008402109
0 6 0.04435162 1 0 591 0.9899552 0.004080157 0.010095589
0 7 0.04449827 1 0 590 0.9882773 0.004404752 0.011791943
0 8 0.04466022 1 0 589 0.9865994 0.004706012 0.013491178
0 9 0.04862678 1 0 588 0.9849215 0.004988167 0.015193306
... 447 more rows
Median duration: 0.5141081
Duration times limited to: 3.595082
Cases: 600 weighted: 600Box 34ehd2p: log S against the residual, with its reference line [book p. 226, Box 8.2.6]
resid_panel <- function(fit, label, file) {
rp <- tda_ple(Surv(Residual, Des) ~ 1, fit$residuals)
rs <- tda_survivor(rp)
rs <- rs[rs$time > 0 & rs$survivor > 0 & rs$survivor < 1, ]
rs$LogSurv <- log(rs$survivor)
p <- tda_ps(rs, file = file, width = 80, height = 40,
xlim = c(0, 4), ylim = c(-3, 0))
p <- tda_pl(p, "plxa", sc = 1, ic = 0)
p <- tda_pl(p, "plya", sc = 1, ic = 0)
p <- tda_pl_lines(tda_pl(p, "plframe"), "time", "LogSurv")
p <- tda_pl_polyline(p, c(0, 3), c(0, -3)) # the standard exponential
tda_pl_text(p, label, at = c(1.5, -0.7))
}
plot(resid_panel(ehd2r, "Exponential Model", "ehd2p.ps"))
The exponential curve sits above the reference at the start, crosses
it around 1.7 and falls away below it after that – the same departure
the transformed survivor function showed in ehh1.
Box 8.2.6’s footnote lists the command files for the other models, and the book’s Figure 8.2.1 puts eight of them side by side. Each is the same three steps on a different fit.
Box 35ehd2p to ehg13p: the pseudoresidual check [book pp. 226-227, Box 8.2.6; Fig. 8.2.1]
resid_fits <- list(
"Exponential Model" = function()
tda_rate(with_x("Surv(TFP, DES)"), d, residuals = TRUE),
"Piecewise Constant Exponential (I)" = function()
pc(with_x("Surv(TFP, DES)"), d, residuals = TRUE),
"Piecewise Constant Exponential (II)" = function()
pc(with_x("Surv(TFP, DES)"), d, model = "exponential_periods2",
tp = "0,24,60", residuals = TRUE),
"Gompertz Model" = function()
fitx("gompertz", on = "xb", residuals = TRUE),
"Weibull Model" = function() fitx("weibull", residuals = TRUE),
"Log-logistic Model" = function() fitx("loglogistic", residuals = TRUE),
"Log-normal Model" = function() fitx("lognormal", residuals = TRUE),
"Sickle Model" = function() fitx("sickle", residuals = TRUE))
# the arrangement of Figure 8.2.1: exponential, Gompertz, Weibull and
# sickle down the left column, the two piecewise constant models,
# log-logistic and log-normal down the right
i <- 0L
panels <- lapply(names(resid_fits), function(nm) {
i <<- i + 1L
resid_panel(resid_fits[[nm]](), nm, sprintf("res%02d.ps", i))
})
r <- tda_combine_ps(split(panels[c(1, 2, 4, 3, 5, 6, 8, 7)], rep(1:4, each = 2)), width = 170, height = 180)
tda_plot_ps(tda_read_ps(file.path(r$dir, "combined.ps")))
The book reads these as the exponential and the Weibull departing most from the reference, the Gompertz, log-logistic and log-normal tracking it closely, and the two piecewise constant models closest of all. They are hints for model selection rather than a goodness-of-fit test: the log-logistic, log-normal and piecewise constant models resemble each other so closely here that the choice between them has to rest on something other than these plots.
9. Semi-parametric Transition Rate Models
Book chapter 9, pp. 229-251. Cox models; command files
ehi1 to ehi13.
ehi1 is the Cox model with the book’s covariates.
tda_coxph() is tda_rate(model = 1) under a
name that says what it is.
Box 36ehi1: a Cox model [book pp. 232-243, Boxes 9.1.2, 9.3.3]
ehi1 <- tda_coxph(with_x("Surv(TFP, DES)"), d)
summary(ehi1)Call: tda_coxph(formula = with_x("Surv(TFP, DES)"), data = d)
Episodes: 600
Model: cox
Converged in 5 iterations
Estimate Std. Error C/Error Pr(>|t|)
EDU 0.06685932 0.02491489 2.6835 0.0072854 **
COHO2 0.41130737 0.11534702 3.5658 0.0003627 ***
COHO3 0.30525142 0.12198338 2.5024 0.0123354 *
LFX -0.00398905 0.00093243 -4.2781 1.885e-05 ***
PNOJ 0.06862675 0.04419255 1.5529 0.1204463
PRES -0.02616777 0.00550043 -4.7574 1.961e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
logLik (starting values): -2584.57
logLik: -2546.776 df: 6 AIC: 5105.551
Null model (exponential, constant rate) logLik: -2514.02ehi1$episodes # the episode table TDA reports first sn org des episodes weighted mean.duration ts.min tf.max excluded
1 1 0 0 142 142 128.17606 0 428 FALSE
2 1 0 1 458 458 49.30349 0 350 FALSEehi2 adds marriage as a time-varying covariate.
time is TDA’s position within the episode, so the
definition has to be evaluated inside TDA rather than computed here;
helpers= writes MDate for it to refer to,
without MDate entering the model.
Box 37ehi2: marriage as a time-varying covariate [book pp. 233-235, Boxes 9.1.3, 9.2.1-9.2.2]
dm <- d
dm$MDate <- ifelse(dm$TMAR <= 0, 10000, dm$TMAR - dm$TStart)
ehi2 <- tda_coxph(stats::reformulate(c(XVARS, "MARR"),
response = "Surv(TFP, DES)"), dm,
helpers = "MDate",
define = list(MARR = "gt(time,MDate)"))
ehi2Call: tda_coxph(formula = stats::reformulate(c(XVARS, "MARR"), response = "Surv(TFP, DES)"),
data = dm, helpers = "MDate", define = list(MARR = "gt(time,MDate)"))
Episodes: 600
Model: cox (TDA 1)
logLik (starting values): -2584.57
logLik: -2546.622
Converged in 5 iterations
Idx SN Org Des MT Variable Coeff Error C/Error Signif
1 1 0 1 A EDU 0.0665 0.0249 2.6647 0.9923
2 1 0 1 A COHO2 0.4056 0.1158 3.5038 0.9995
3 1 0 1 A COHO3 0.2901 0.1249 2.3225 0.9798
4 1 0 1 A LFX -0.0042 0.0010 -4.1789 1.0000
5 1 0 1 A PNOJ 0.0674 0.0443 1.5220 0.8720
6 1 0 1 A PRES -0.0263 0.0055 -4.7804 1.0000
7 1 0 1 A MARR 0.0660 0.1192 0.5541 0.4205ehi2$episodes # the episode table TDA reports first sn org des episodes weighted mean.duration ts.min tf.max excluded
1 1 0 0 142 142 128.17606 0 428 FALSE
2 1 0 1 458 458 49.30349 0 350 FALSEehi3, ehi4 and
ehi5 are Cox models on the split episodes of
ehf – covariates alone, then with marriage, then
distinguishing married men.
Box 38ehi3, ehi4, ehi5: Cox models on the split episodes [book pp. 237-239, Boxes 9.2.3-9.2.6]
ehi3 <- tda_coxph(stats::reformulate(XVARS,
response = "Surv(ts, tf, des)"), sp)
ehi3Call: tda_coxph(formula = stats::reformulate(XVARS, response = "Surv(ts, tf, des)"),
data = sp)
Episodes: 761
Model: cox (TDA 1)
logLik (starting values): -2584.57
logLik: -2546.776
Converged in 5 iterations
Idx SN Org Des MT Variable Coeff Error C/Error Signif
1 1 0 1 A EDU 0.0669 0.0249 2.6835 0.9927
2 1 0 1 A COHO2 0.4113 0.1153 3.5658 0.9996
3 1 0 1 A COHO3 0.3053 0.1220 2.5024 0.9877
4 1 0 1 A LFX -0.0040 0.0009 -4.2781 1.0000
5 1 0 1 A PNOJ 0.0686 0.0442 1.5529 0.8796
6 1 0 1 A PRES -0.0262 0.0055 -4.7574 1.0000ehi3$episodes # all three are fitted to these same split episodes sn org des episodes weighted mean.duration ts.min tf.max excluded
1 1 0 0 303 303 76.35313 0 428 FALSE
2 1 0 1 458 458 38.53057 0 350 FALSEtda_coxph(stats::reformulate(c(XVARS, "MARR"),
response = "Surv(ts, tf, des)"), sp)Call: tda_coxph(formula = stats::reformulate(c(XVARS, "MARR"), response = "Surv(ts, tf, des)"),
data = sp)
Episodes: 761
Model: cox (TDA 1)
logLik (starting values): -2584.57
logLik: -2546.622
Converged in 5 iterations
Idx SN Org Des MT Variable Coeff Error C/Error Signif
1 1 0 1 A EDU 0.0665 0.0249 2.6647 0.9923
2 1 0 1 A COHO2 0.4056 0.1158 3.5038 0.9995
3 1 0 1 A COHO3 0.2901 0.1249 2.3225 0.9798
4 1 0 1 A LFX -0.0042 0.0010 -4.1789 1.0000
5 1 0 1 A PNOJ 0.0674 0.0443 1.5220 0.8720
6 1 0 1 A PRES -0.0263 0.0055 -4.7804 1.0000
7 1 0 1 A MARR 0.0660 0.1192 0.5541 0.4205tda_coxph(stats::reformulate(c(XVARS, "MARR", "MarrMen"),
response = "Surv(ts, tf, des)"), sp)Call: tda_coxph(formula = stats::reformulate(c(XVARS, "MARR", "MarrMen"),
response = "Surv(ts, tf, des)"), data = sp)
Episodes: 761
Model: cox (TDA 1)
logLik (starting values): -2584.57
logLik: -2528.446
Converged in 5 iterations
Idx SN Org Des MT Variable Coeff Error C/Error Signif
1 1 0 1 A EDU 0.0792 0.0241 3.2867 0.9990
2 1 0 1 A COHO2 0.3767 0.1158 3.2537 0.9989
3 1 0 1 A COHO3 0.2928 0.1243 2.3556 0.9815
4 1 0 1 A LFX -0.0045 0.0010 -4.4408 1.0000
5 1 0 1 A PNOJ 0.1012 0.0449 2.2524 0.9757
6 1 0 1 A PRES -0.0252 0.0054 -4.6858 1.0000
7 1 0 1 A MARR 0.5046 0.1303 3.8733 0.9999
8 1 0 1 A MarrMen -0.8608 0.1420 -6.0624 1.0000ehi7 and ehi8 ask the proportionality question twice: once as an interaction with log time, once by letting the effect vary between periods. Either way, a proportional effect is one that does not move.
Box 39ehi7, ehi8: proportional hazards [book pp. 243-245, Boxes 9.3.2, 9.3.4-9.3.5]
dw <- d
dw$Women <- as.integer(dw$SEX == 2)
ehi7 <- tda_coxph(stats::reformulate(c(XVARS, "Women", "WTest"),
response = "Surv(TFP, DES)"), dw,
define = list(WTest = "Women * (log(time) - 4.73)"))
ehi7Call: tda_coxph(formula = stats::reformulate(c(XVARS, "Women", "WTest"),
response = "Surv(TFP, DES)"), data = dw, define = list(WTest = "Women * (log(time) - 4.73)"))
Episodes: 600
Model: cox (TDA 1)
logLik (starting values): -2584.57
logLik: -2536.272
Converged in 5 iterations
Idx SN Org Des MT Variable Coeff Error C/Error Signif
1 1 0 1 A EDU 0.0803 0.0247 3.2566 0.9989
2 1 0 1 A COHO2 0.3937 0.1155 3.4097 0.9993
3 1 0 1 A COHO3 0.2716 0.1222 2.2231 0.9738
4 1 0 1 A LFX -0.0041 0.0009 -4.3657 1.0000
5 1 0 1 A PNOJ 0.0888 0.0447 1.9872 0.9531
6 1 0 1 A PRES -0.0259 0.0054 -4.7654 1.0000
7 1 0 1 A Women 0.7111 0.1607 4.4244 1.0000
8 1 0 1 A WTest 0.2587 0.0987 2.6211 0.9912ehi7$episodes sn org des episodes weighted mean.duration ts.min tf.max excluded
1 1 0 0 142 142 128.17606 0 428 FALSE
2 1 0 1 458 458 49.30349 0 350 FALSEehi8 <- tda_coxph(stats::reformulate(c(XVARS, "Women"),
response = "Surv(TFP, DES)"), dw,
options = list(tp = "0 (12) 96"))
summary(ehi8) # the global goodness-of-fit test comes with the summaryCall: tda_coxph(formula = stats::reformulate(c(XVARS, "Women"), response = "Surv(TFP, DES)"),
data = dw, options = list(tp = "0 (12) 96"))
Episodes: 600
Model: cox
Converged in 5 iterations
Estimate Std. Error C/Error Pr(>|t|)
EDU 0.07632115 0.02466457 3.0944 0.0019724 **
COHO2 0.38596002 0.11536101 3.3457 0.0008208 ***
COHO3 0.29586568 0.12199590 2.4252 0.0152995 *
LFX -0.00408664 0.00093302 -4.3800 1.187e-05 ***
PNOJ 0.09030560 0.04482867 2.0145 0.0439612 *
PRES -0.02492391 0.00542188 -4.5969 4.288e-06 ***
Women 0.36988945 0.09763713 3.7884 0.0001516 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
logLik (starting values): -2584.57
logLik: -2539.679 df: 7 AIC: 5093.358
Null model (exponential, constant rate) logLik: -2514.02
Global goodness-of-fit:
sn org des statistic df significance
1 0 1 62.92011 56 0.7552702ehi6 asks it graphically instead:
log(-log(S)) per sex, parallel if the hazards are
proportional. The axes, tick spacing and the two labels are
ehi6.cf’s own.
Box 40ehi6: the graphical check [book pp. 241-242, Box 9.3.1; Fig. 9.3.1]
lls <- tda_survivor(pl2, conf.int = NULL)
lls <- lls[lls$survivor > 0 & lls$survivor < 1, ]
lls$LLS <- log(-log(lls$survivor))
p <- tda_ps(lls, width = 80, height = 40, xlim = c(0, 240), ylim = c(-5, 1))
p <- tda_pl(p, "plxa", sc = 24, ic = 0)
p <- tda_pl(p, "plya", sc = 1, ic = 0)
p <- tda_pl_lines(p, "time", "LLS", by = "group", lty = 1) # ehi6.cf sets no lt=
p <- tda_pl_text(p, "Women", at = c(192, 0.6))
plot(tda_pl_text(p, "Men", at = c(216, 0.1)))
ehi9 gives each sex its baseline rate instead of assuming one.
Box 41ehi9: a stratified Cox model [book pp. 246-247, Boxes 9.3.6-9.3.7]
ehi9 <- tda_coxph(stats::reformulate(c(XVARS, "strata(SEX)"),
response = "Surv(TFP, DES)"), d)
ehi9Call: tda_coxph(formula = stats::reformulate(c(XVARS, "strata(SEX)"),
response = "Surv(TFP, DES)"), data = d)
Episodes: 600
Model: cox (TDA 1)
logLik (starting values): -2259.567
logLik: -2224.632
Converged in 5 iterations
Idx SN Org Des MT Variable Coeff Error C/Error Signif
1 1 0 1 A EDU 0.0793 0.0248 3.2013 0.9986
2 1 0 1 A COHO2 0.3975 0.1156 3.4391 0.9994
3 1 0 1 A COHO3 0.2801 0.1226 2.2848 0.9777
4 1 0 1 A LFX -0.0040 0.0009 -4.3158 1.0000
5 1 0 1 A PNOJ 0.0889 0.0448 1.9856 0.9529
6 1 0 1 A PRES -0.0258 0.0054 -4.7261 1.0000ehi9$episodes # one block per stratum, as TDA reports them sn org des episodes weighted mean.duration ts.min tf.max excluded
1 1 0 0 142 142 128.17606 0 428 FALSE
2 1 0 1 458 458 49.30349 0 350 FALSEehi10 evaluates the baseline rate and its cumulation, ehi11 draws it, ehi12 smooths it with a spline and ehi13 draws the smoothed rate.
Box 42ehi10: the baseline rate, two cohorts [book pp. 249-250, Boxes 9.4.1-9.4.2]
# ehi10.cf asks for two constellations, cohort 1 and cohort 3, both
# written to ehi10.bl; they come back as sub-tables of $rates, told
# apart by ID
ehi10 <- tda_coxph(with_x("Surv(TFP, DES)"), d,
prate = list(
list(EDU = 13, COHO2 = 0, COHO3 = 0,
LFX = 5, PNOJ = 1, PRES = 30),
list(EDU = 13, COHO2 = 0, COHO3 = 1,
LFX = 5, PNOJ = 1, PRES = 30)))
ehi10Call: tda_coxph(formula = with_x("Surv(TFP, DES)"), data = d, prate = list(list(EDU = 13,
COHO2 = 0, COHO3 = 0, LFX = 5, PNOJ = 1, PRES = 30), list(EDU = 13,
COHO2 = 0, COHO3 = 1, LFX = 5, PNOJ = 1, PRES = 30)))
Episodes: 600
Model: cox (TDA 1)
logLik (starting values): -2584.57
logLik: -2546.776
Converged in 5 iterations
Idx SN Org Des MT Variable Coeff Error C/Error Signif
1 1 0 1 A EDU 0.0669 0.0249 2.6835 0.9927
2 1 0 1 A COHO2 0.4113 0.1153 3.5658 0.9996
3 1 0 1 A COHO3 0.3053 0.1220 2.5024 0.9877
4 1 0 1 A LFX -0.0040 0.0009 -4.2781 1.0000
5 1 0 1 A PNOJ 0.0686 0.0442 1.5529 0.8796
6 1 0 1 A PRES -0.0262 0.0055 -4.7574 1.0000ehi10$episodes # the episodes behind the baseline sn org des episodes weighted mean.duration ts.min tf.max excluded
1 1 0 0 142 142 128.17606 0 428 FALSE
2 1 0 1 458 458 49.30349 0 350 FALSEbl <- ehi10$rates
cols <- c("ID", "Time", "RiskSet", "CumRate", "BaselineRate")
head(bl[, cols], 10) ID Time RiskSet CumRate BaselineRate
1 0 0 600 0.000000000 0.000000000
2 0 2 600 0.004324803 0.001893449
3 0 3 597 0.015201628 0.009524000
4 0 4 590 0.034991229 0.017328233
5 0 5 581 0.041721469 0.005893155
6 0 6 577 0.064307187 0.019776577
7 0 7 567 0.085179004 0.018275846
8 0 8 557 0.099372708 0.012428336
9 0 9 548 0.116206159 0.014739759
10 0 10 540 0.135804437 0.017160707tail(bl[bl$ID == 0, cols], 8) ID Time RiskSet CumRate BaselineRate
124 0 215 39 2.811890 0.0068633098
125 0 220 38 2.851719 0.0069749534
126 0 275 26 2.909198 0.0009150944
127 0 293 20 2.986824 0.0037761719
128 0 312 16 3.076748 0.0041441830
129 0 326 14 3.176270 0.0062245309
130 0 332 11 3.298348 0.0178158501
131 0 350 9 3.443908 0.0070808389ehi11.cf draws the two cumulated rates against time (the
book’s Figure 9.4.1).
Box 43ehi11: the two cohorts’ cumulated rates [book p. 250, Fig. 9.4.1]
p <- tda_ps(bl, width = 80, height = 40, xlim = c(0, 300), ylim = c(0, 5))
p <- tda_pl(p, "plxa", sc = 60, ic = 5)
p <- tda_pl(p, "plya", sc = 1, ic = 0)
p <- tda_pl_lines(p, "Time", "CumRate", by = "ID", select = "0", lty = 1)
p <- tda_pl_lines(p, "Time", "CumRate", by = "ID", select = "1", lty = 5)
p <- tda_pl_text(p, "COHO1", at = c(240, 2.3))
plot(tda_pl_text(p, "COHO3", at = c(240, 4.2)))
The rate itself jumps at every event time, so the book does not
smooth it directly. ehi12 smooths the cumulative
rate – column 7 of ehi10.bl, the one for cohort 1 – with a
cubic spline over a half-month grid, and ehi13 plots that
spline’s first derivative, which is the rate.
Box 44ehi12: the cumulated rate smoothed [book p. 251]
# ehi12.cf selects the first sub-table, cohort 1
b1 <- bl[bl$ID == 0, ]
sm <- tda_spl(b1$Time, b1$CumRate, sig = 0.2, deg = 3,
rx = seq(0, 300, by = 0.5))$table
head(sm[, c("x", "fitted", "d1")], 4) x fitted d1
1 0.0 -0.026485261 0.01191835
2 0.5 -0.020365553 0.01255804
3 1.0 -0.013929678 0.01318301
4 1.5 -0.007184992 0.01379328ehi13.cf draws the smoother’s first derivative, which is
the rate itself (the book’s Figure 9.4.2).
Box 45ehi13: the smoothed baseline rate [book p. 251, Fig. 9.4.2]
10. Models with a Mixture Distribution
Book chapter 10, pp. 261-271. A gamma-distributed error term on the exponential and the Weibull; no command file ships for either, so these two fits are checked against the book’s boxes rather than against TDA’s reference output.
The book’s chapter 10 adds unobserved heterogeneity:
mix = 1 is a gamma mixture on the exponential model (Box
10.2.1) and on the Weibull (Box 10.2.2); D Constant is the
log of the mixture variance.
Box 46Boxes 10.2.1, 10.2.2: gamma mixtures [book pp. 268-271, Boxes 10.2.1-10.2.2]
tda_rate(with_x("Surv(TFP, DES)"), d, options = list(mix = 1))Call: tda_rate(formula = with_x("Surv(TFP, DES)"), data = d, options = list(mix = 1))
Episodes: 600
Model: exponential (TDA 2)
logLik (starting values): -2521.785
logLik: -2442.128
Converged in 6 iterations
Idx SN Org Des MT Variable Coeff Error C/Error Signif
1 1 0 1 A Constant -3.9869 0.3430 -11.6219 1.0000
2 1 0 1 A EDU 0.0873 0.0310 2.8108 0.9951
3 1 0 1 A COHO2 0.5707 0.1445 3.9503 0.9999
4 1 0 1 A COHO3 0.4163 0.1512 2.7529 0.9941
5 1 0 1 A LFX -0.0046 0.0011 -4.2132 1.0000
6 1 0 1 A PNOJ 0.0850 0.0545 1.5594 0.8811
7 1 0 1 A PRES -0.0321 0.0066 -4.8423 1.0000
8 1 0 1 D Constant -0.8019 0.1837 -4.3659 1.0000tda_rate(with_x("Surv(TFP, DES)"), d, model = "weibull",
options = list(mix = 1))Call: tda_rate(formula = with_x("Surv(TFP, DES)"), data = d, model = "weibull",
options = list(mix = 1))
Episodes: 600
Model: weibull (TDA 7)
logLik (starting values): -2521.785
logLik: -2410.713
Converged in 6 iterations
Idx SN Org Des MT Variable Coeff Error C/Error Signif
1 1 0 1 A Constant -3.4926 0.2953 -11.8272 1.0000
2 1 0 1 A EDU 0.0720 0.0266 2.7106 0.9933
3 1 0 1 A COHO2 0.4641 0.1205 3.8505 0.9999
4 1 0 1 A COHO3 0.4014 0.1284 3.1274 0.9982
5 1 0 1 A LFX -0.0041 0.0008 -4.8731 1.0000
6 1 0 1 A PNOJ 0.1144 0.0470 2.4345 0.9851
7 1 0 1 A PRES -0.0260 0.0053 -4.8769 1.0000
8 1 0 1 B Constant 0.6693 0.0880 7.6057 1.0000
9 1 0 1 D Constant 0.7833 0.1716 4.5648 1.0000Finally, the product-limit survivor function on its own, the picture
the whole exercise starts from. This one answers to no command file: the
book has no such figure, and ple5p.cf – which the name
would suggest – draws something else, the two groups with their
confidence bands, as ehc8 already does above.
Box 47The product-limit survivor function (no book box)
p <- tda_ps(tda_survivor(pl), width = 80, height = 40, xlim = c(0, 300), ylim = c(0, 1))
p <- tda_pl(p, "plxa", sc = 60, ic = 5)
p <- tda_pl(p, "plya", sc = 1, ic = 10)
plot(tda_pl_lines(p, "time", "survivor"))