How can I use survival analysis with multi-state model to look at a progressive events (gentrification stages)?

Question

I'm working examining the effect of a particular anti-gentrification policy in Berlin. I want to use a survival analysis to see if the amount of time it takes any area to progress to the next stage of gentrification out of five stages (ungentrified through gentrified), but a normal regression doesn't work here because some areas might stop gentrifying altogether, and thus have 0 or NA for the time to the next stage. I've tried using R's survival and mstate packages but haven't had much luck finding relevant literature for this sort of analysis, and I'm wondering if anyone has any suggestions on how best to go about this.

For some further context, here's some of the dataset I'm using, where gcode is the final status in 2019, and tg0:tg5 is the time it takes to get from stage 0 to 1, 1 to 2, and so on. As you can see, the 0.000 places represent censored data where an area never reached the next stage.

> head(data, n=20)
    RAUMID gcode   tg0    tg1   tg2   tg3   tg4   tg5
1  1011101     4 2.001  1.999 2.001 4.000 0.000    NA
2  1011102     5 0.999  4.000 4.000 2.001 1.999 5.002
3  1011103     2 0.999 12.000 0.000 0.000 0.000    NA
4  1011104     4 0.999  6.001 1.999 4.000 0.000    NA
5  1011105     4 0.999  4.000 2.001 9.999 0.000    NA
6  1011201     4 0.999  6.001 1.999 6.001 0.000    NA
7  1011202     4 0.999  2.001 1.999 6.001 0.000    NA
8  1011203     4 0.999  4.000 2.001 4.000 0.000    NA
9  1011204     4 0.999  6.001 1.999 2.001 0.000    NA
10 1011301     2 0.999 12.000 0.000 0.000 0.000    NA
11 1011302     1 0.999  0.000 0.000 0.000 0.000    NA
12 1011303     5 0.999  2.001 1.999 2.001 4.000 7.001
13 1011304     2 0.000  0.999 0.000 0.000 0.000    NA
14 1011305     4 0.999  8.000 6.001 1.999 0.000    NA
15 1011306     2 0.000  0.999 0.000 0.000 0.000    NA
16 1011401     2 2.001  1.999 0.000 0.000 0.000    NA
17 1011402     1 2.001  0.000 0.000 0.000 0.000    NA
18 1022101     4 0.000  0.999 4.000 4.000 0.000    NA
19 1022102     4 0.000  0.999 4.000 4.000 0.000    NA
20 1022103     2 0.000 11.001 0.000 0.000 0.000    NA

Creating dummy variables for whether an area has been censored or whether it was treated during the period when it transitioned from one stage to another has been no problem, but actually using survival analysis for this type of problem has been really tricky.

Edit: for some additional context, here is a sample of data from which I calculated tg0 through tg5:

To provide some further context, stage 0 means not (yet) gentrifiable, and is calculated based on an area being above the 70th percentile in average sale prices for residential properties. While I went back and forth on whether these areas should simply be eliminated from the dataset entirely before continued categorization, I decided to include them as, theoretically, an area could be initially ungentrifiable but then decrease in value sufficiently over the course of 10 years to be gentrifiable.

> head(data, n=20)
    RAUMID   gcode0yr   gcode1yr   gcode2yr   gcode3yr   gcode4yr   gcode5yr
1  1011101 2002-12-31 2004-12-31 2006-12-31 2008-12-31 2012-12-31       <NA>
2  1011102 2001-12-31 2002-12-31 2006-12-31 2010-12-31 2012-12-31 2014-12-31
3  1011103 2001-12-31 2002-12-31 2014-12-31       <NA>       <NA>       <NA>
4  1011104 2001-12-31 2002-12-31 2008-12-31 2010-12-31 2014-12-31       <NA>
5  1011105 2001-12-31 2002-12-31 2006-12-31 2008-12-31 2018-12-31       <NA>
6  1011201 2001-12-31 2002-12-31 2008-12-31 2010-12-31 2016-12-31       <NA>
7  1011202 2001-12-31 2002-12-31 2004-12-31 2006-12-31 2012-12-31       <NA>
8  1011203 2001-12-31 2002-12-31 2006-12-31 2008-12-31 2012-12-31       <NA>
9  1011204 2001-12-31 2002-12-31 2008-12-31 2010-12-31 2012-12-31       <NA>
10 1011301 2001-12-31 2002-12-31 2014-12-31       <NA>       <NA>       <NA>
11 1011302 2001-12-31 2002-12-31       <NA>       <NA>       <NA>       <NA>
12 1011303 2001-12-31 2002-12-31 2004-12-31 2006-12-31 2008-12-31 2012-12-31
13 1011304       <NA> 2001-12-31 2002-12-31       <NA>       <NA>       <NA>
14 1011305 2001-12-31 2002-12-31 2010-12-31 2016-12-31 2018-12-31       <NA>
15 1011306       <NA> 2001-12-31 2002-12-31       <NA>       <NA>       <NA>
16 1011401 2002-12-31 2004-12-31 2006-12-31       <NA>       <NA>       <NA>
17 1011402 2002-12-31 2004-12-31       <NA>       <NA>       <NA>       <NA>
18 1022101       <NA> 2001-12-31 2002-12-31 2006-12-31 2010-12-31       <NA>
19 1022102       <NA> 2001-12-31 2002-12-31 2006-12-31 2010-12-31       <NA>
20 1022103       <NA> 2001-12-31 2012-12-31       <NA>       <NA>       <NA>

Another edit: Here is the outcome of my data wrangling:

 dataL %>% group_by(RAUMID) %>% head(n=10)
# A tibble: 10 x 6
# Groups:   RAUMID [10]
    RAUMID startTime endTime event treatment_starttime treatedduring
     <dbl>     <dbl>   <dbl> <fct>               <dbl>         <dbl>
 1 1011101      4       6.00 tg2                  NA               0
 2 1011102      2.00    6.00 tg2                  NA               0
 3 1011103      2.00   14.0  tg2                  18.0             0
 4 1011104      2.00    8    tg2                  NA               0
 5 1011105      2.00    6.00 tg2                  NA               0
 6 1011201      2.00    8    tg2                  NA               0
 7 1011202      2.00    4    tg2                  NA               0
 8 1011203      2.00    6.00 tg2                  NA               0
 9 1011204      2.00    8    tg2                  NA               0
10 1011301      2.00   14.0  tg2                  NA               0
> dataL %>% group_by(RAUMID) %>% filter(treatedduring == 1) %>%  head(n=10)
# A tibble: 10 x 6
# Groups:   RAUMID [10]
    RAUMID startTime endTime event treatment_starttime treatedduring
     <dbl>     <dbl>   <dbl> <fct>               <dbl>         <dbl>
 1 2020204     0.999    6.00 tg2                 2.18              1
 2 2020206     0.999    2.00 tg2                -5.43              1
 3 2030401     0.999    4    tg2                -5.43              1
 4 2050802     4        6.00 tg2                -1.71              1
 5 3040614     2.00    19.0  cens               -0.791             1
 6 3050925     2.00    19.0  cens               16.9               1
 7 3061126     0.999   19.0  cens               -1.74              1
 8 3061131     4       19.0  cens               -3.78              1
 9 3061227     0.999   19.0  cens               -0.153             1
10 3061332     2.00    19.0  cens               13.5               1

Answer 1

The multi-state vignette for the R survival package sets out how to approach this type of problem; your situation with a set of sequential, uni-directional possible transitions is illustrated in the top right panel of Figure 1 of the vignette. So in principle your situation can be handled with a survival model.

The main problems here are: (1) formatting the data correctly, (2) formulating the model appropriately and (3) geting a useful result even if you do all that. Item (3) can't really be predicted in advance. Here's what you need to watch out for in Items (1) and (2).

The main task for Item (1) is to convert these wide-form data into long-form data, with each row representing a particular defined time period for a single RAUMID. Each row then specifies the RAUMID, a start time, a stop time, any covariate values in place during the period (in particular, was the anti-gentrification policy in place) and whether there was a change of state (event) at the end of the period.

That (startTime,stopTime,event) formulation requires two additional steps starting from your data. First, the times you currently show are for the durations within each state. You probably should convert those all to time elapsed since study start (taken to be time = 0).

Second, you need to specify a categorical variable that represents an event, a change in state in your situation. With the survival package the reference level of that variable needs to represent censoring, or no change in state at stopTime. As there are no transitions into tg0 you don't need to include a level for that, but you would need to specify categories for tg1, tg2, tg3, tg4, and tg5. For each (startTime, stopTime, event), if there is a change of state at stopTime you set the event variable to the corresponding new category. That format also allows for a change in covariate values over time; if the time period in question didn't end with an event change, the event category is listed as the reference level representing censoring. If I'm interpreting your data display correctly, the first row of your data would be changed into multiple rows as follows:

RAUMID	startTime	stopTime	event	policy
1011101	0	2.001	tg1	p?
1011101	2.001	4.000	tg2	p?
1011101	4.000	6.001	tg3	p?
1011101	6.001	10.001	tg4	p?
1011101	10.001	t?	cens	p?

The t? in the last stopTime should be replaced by the time elapsed between time = 0 and the last follow-up, when the RAUM was still in state tg4. The policy value, shown as p? above, represents the policy in place during that time interval. You append whatever other covariate values held during each time period within the RAUM as additional columns. If a new policy comes into play at time $T_p$ after time = 0, set a stopTime equal to $T_p$ and a censoring event for all affected RAUMID (unless there was also an event change at that time) and start a new data row with startTime of $T_p$.

Your sample data indicate a couple of things that you'll need to clarify. First, some RAUMID show values for tg5, which I interpret to mean time spent in tg5. If there is no transition possible out of tg5 it's not clear what that's intended to represent. Also, some RAUMID show values of 0 for the time in state tg0. Presumably those are cases in which the first observed state was tg1. If that's the case, you would specify the initial state with an istate argument ; see page 8 of the multi-state vignette.

The main survival package vignette provides guidance on tools available to help with this type of data formatting, and in particular the survcheck() function that helps identify the errors that are almost inevitable in practice. Once the data are in that format, the multi-state vignette shows how to generate curves showing things like the probability in each state as a function of time

For Item (2), you presumably want to see if policy changes affected the rate of transitions. You then need to apply your knowledge of the subject matter, in particular: do you expect the change in policy to affect all transitions similarly, or might it have different effects on different transitions? Section 3.4.2 of the main survival vignette illustrates how you can allow for either common or different effects of covariates on each transition.

As the data are only collected annually, you will probably be better off approaching this as a discrete-time survival analysis. That can be done as an appropriately formatted logistic regression. John Willett has a website with several references to publications with his colleagues, in particular Judith Singer, showing how this can be done with logistic regression. Paul Allison wrote a characteristically clear review in 1982 discussing the principles and describing both logistic regression and complementary log-log modeling approaches. Tutz and Schmid have a more recent book-length treatment with R code, including how to use tools provided by the discSurv package.

One possibility would be to format your data with one row per RAUMID per year, including (along with other covariates) the starting tg state as a categorical predictor, whether the policy was in effect that year, and a 0/1 outcome variable for no transition versus transition to next tg state that year. Logistic regression then would provide coefficients representing the effects of starting tg state and of policy on the log-odds of transition to the next state. You could include an interaction term between tg and policy in the model to see whether the policy has different associations with transitions depending on the tg state.

Finally, from comments it seems that your assignments to tg states is based on a categorization process done every 2 years. That has two potential implications. For one, if you only do the categorization every 2 years, you might need to re-cast your data format into 2-year bins. For another, you might consider whether you would be better off primarily modeling a continuous measure of the gentrification process, then perhaps translating the results with the continuous measure into your categorization.