Plot the estimated hazard curve at mean-1*SD mean, and mean+1*SD after fitting survival analysis using Stata.

[smith Jason]

Hi, everybody,
I have a dataset like this,
* Example generated by -dataex.
clear
input int id int time byte injury byte safescore
1 1 1 9
2 1 1 9
3 1 0 6
3 2 0 6
3 3 1 6
4 1 0 6
4 2 0 6
4 3 0 6
5 1 0 6
5 2 1 7
6 1 0 6
6 2 0 6
7 1 0 6
7 2 0 8
8 1 0 6
8 2 0 6
8 3 1 7
9 1 0 6
9 2 1 8
10 1 0 6
10 2 1 7
end

For ease of interpretation, I centered the safescore at the mean and created the resulting variable c_safescore.
Meanwhile, I created the corresponding mean+1*SD and mean-1*SD variables as follows,

egen safescore_mean=mean(safescore)
egen safescore_sd=sd(safescore)
gen c_safescore=safescore-safescore_mean
gen c_safescore_upper=c_safescore+1*safescore_sd
gen c_safescore_lower=c_safescore-1*safescore_sd

I don't know if the following procedures are incorrect to achieve my goal:

To achieve what I want, I fitted a simple discrete-time survival model with logit function in Stata.
logit injury c_safescore,nocons
logit injury c_safescore_upper,nocons
logit injury c_safescore_lower,nocons

Now, I want to graph the fitted hazard curve at mean of safescore, mean+SD, and mean-SD in Stata in a plane.

Can someone help me to do this in Stata?

My own solution is below,

predict c_safescore_m
line c_safescore_m time, sort
||
line c_safescore_lower time, sort
||
line c_safescore_upper time, sort

It seems that something is incorrect.

[smith Jason]

My question is I don't know if I should use c_safescore as mean value of safescore to plot or use mean(safescore).
I am confused about it and don't know how to graph it.

[smith Jason]

anybody?

[Joseph Coveney]

To achieve what I want, I fitted a simple discrete-time survival model with logit function in Stata. . . . Now, I want to graph the fitted hazard curve at mean of safescore, mean+SD, and mean-SD

You might be looking for something more like the following.

bysort id (time): generate byte last = _n == _N
quietly summarize safescore if last
local mmsd = r(mean) - r(sd)
local mpsd = r(mean) + r(sd)
local m = r(mean)

logit injury i.time c.safescore, or nolog
margins time, at(safescore = (`mmsd' `m' `mpsd'))
marginsplot , xdimension(time) title("") ///
    ylabel( , format(%03.2f) angle(horizontal) nogrid) ///
    noci legend(off)  plotopts(msymbol(none))

[smith Jason]

Thank you!

[smith Jason]

You might be looking for something more like the following.

bysort id (time): generate byte last = _n == _N
quietly summarize safescore if last
local mmsd = r(mean) - r(sd)
local mpsd = r(mean) + r(sd)
local m = r(mean)

logit injury i.time c.safescore, or nolog
margins time, at(safescore = (`mmsd' `m' `mpsd'))
marginsplot , xdimension(time) title("") ///
ylabel( , format(%03.2f) angle(horizontal) nogrid) ///
noci legend(off) plotopts(msymbol(none))

Why don't you used centered safescore variable in the logit model?
Also, I don't understand why don't you use the overall mean of safescore variable?
Thank you!

[Joseph Coveney]

Why don't you used centered safescore variable in the logit model?
Also, I don't understand why don't you use the overall mean of safescore variable?

As to not using the overall average, the individual values are not independent within ID, and there are an unequal number of observations between IDs. So, some IDs would contribute their presumably correlated values of safescore more than others, making the pooled-value mean representative of Idontknow. In order to help assure that independent observations contribute to the mean, I chose each ID's last observation inasmuch as it is the closest in time to the injury or its most recently observed avoidance.

If you want to get fancy about using some kind of grand mean, then I suppose that you can do something like the following in recognition that the observations are not independent.

mixed safescore i.tim || id: , reml dfmethod(kroger) nolrtest nolog
 
tempname mean SD
estat sd
scalar define `SD' = r(b)[1, "Residual:sd(e)"]
 
margins , atmeans df(`e(ddf_m)')
scalar define `mean' = r(table)["b", "_cons"]
 
local mmsd = `mean' - `SD'
local mpsd = `mean' + `SD'
local m =  `mean'

You'd use those local macro variables in place of the ones from yesterday.

But, the purpose is to display the relationship between variations in safescore and variations in hazard, and for that you can just as easily choose three values of safescore arbitrarily. The values of safescore in your dataset range only from six to nine; you could have just chosen the minimum, maximum and a midway value to illustrate the point you want to make on the hazard plot.

As to centering, I'm not sure of what your purpose is in wanting to center the values. Or what would you use to center them (grand mean of whatever flavor?, each individual ID’s mean?).

[smith Jason]

Below is the updated dataset for demonstration.
* Example generated by -dataex.
clear
input int id int time byte injury byte instruction byte safescore
1 1 1 0 8
2 1 1 0 7
3 1 0 0 6
3 2 0 0 5
3 3 0 0 4
3 4 0 0 3
3 5 1 0 2
4 1 0 0 1
4 2 0 0 1
4 3 0 0 2
4 4 0 0 3
4 5 0 0 4
5 1 0 0 5
5 2 0 0 6
5 3 0 0 7
5 4 0 0 7
5 5 1 0 5
6 1 0 0 4
6 2 0 0 3
6 3 0 0 2
6 4 0 0 6
6 5 1 0 5
7 1 0 0 4
7 2 0 0 3
7 3 0 0 2
7 4 0 0 4
7 5 1 0 3
8 1 0 0 7
8 2 0 0 8
8 3 0 0 4
8 4 0 0 5
8 5 1 0 6
9 1 0 0 3
9 2 0 0 2
9 3 0 0 6
9 4 0 0 3
9 5 0 0 9
9 6 1 0 1
10 1 0 0 2
10 2 0 0 5
10 3 0 0 6
10 4 0 0 7
10 5 0 0 8
10 6 1 0 3
11 1 0 0 2
11 2 0 0 1
11 3 0 0 4
11 4 0 0 3
11 5 0 0 7
11 6 1 0 2
12 1 0 0 7
12 2 0 0 5
12 3 0 0 4
12 4 0 0 3
12 5 1 0 2
13 1 0 0 5
13 2 0 0 6
13 3 0 0 7
13 4 0 0 2
13 5 0 0 4
13 6 0 0 6
13 7 0 0 7
14 1 0 0 7
14 2 0 0 4
14 3 0 0 3
14 4 0 0 2
14 5 0 0 4
14 6 0 0 5
15 1 0 1 9
15 2 0 1 8
15 3 0 1 7
15 4 0 1 6
15 5 1 1 5
16 1 0 1 4
16 2 0 1 3
16 3 0 1 2
16 4 0 1 6
16 5 1 1 8
17 1 0 1 9
17 2 0 1 3
17 3 0 1 1
17 4 0 1 8
17 5 1 1 4
18 1 0 1 3
18 2 0 1 2
18 3 0 1 5
18 4 0 1 5
18 5 1 1 2
19 1 0 1 1
19 2 0 1 3
19 3 0 1 8
19 4 0 1 6
19 5 1 1 3
20 1 0 1 5
20 2 0 1 4
20 3 0 1 3
20 4 0 1 2
20 5 1 1 1
21 1 0 1 7
21 2 0 1 5
21 3 0 1 6
21 4 0 1 7
21 5 0 1 9
21 6 1 1 1
22 1 0 1 2
22 2 0 1 3
22 3 0 1 4
22 4 0 1 5
22 5 0 1 8
22 6 1 1 2
23 1 0 1 9
23 2 0 1 3
23 3 0 1 5
23 4 0 1 2
23 5 0 1 4
24 1 0 1 8
24 2 0 1 4
24 3 0 1 3
24 4 0 1 7
24 5 1 1 1
25 1 0 1 2
25 2 0 1 6
25 3 0 1 4
25 4 0 1 3
25 5 0 1 5
25 6 0 1 2
26 1 0 1 2
26 2 0 1 5
26 3 0 1 6
26 4 0 1 7
26 5 0 1 3
26 6 0 1 2
27 1 0 1 4
27 2 0 1 3
27 3 0 1 6
27 4 0 1 9
27 5 0 1 6
27 6 1 1 5
28 1 0 1 4
28 2 0 1 3
28 3 0 1 2
28 4 0 1 1
28 5 1 1 9
29 1 0 1 1
29 2 0 1 8
29 3 0 1 2
29 4 0 1 3
29 5 1 1 4
30 1 0 1 5
30 2 0 1 6
30 3 0 1 7
30 4 0 1 8
30 5 1 1 8
31 1 0 1 2
31 2 0 1 5
31 3 1 1 6
32 1 0 1 6
32 2 0 1 7
32 3 0 1 3
32 4 0 1 5
32 5 0 1 6
32 6 0 1 4
32 7 1 1 2
33 1 0 1 5
33 2 0 1 6
33 3 0 1 5
33 4 0 1 7
33 5 0 1 7
33 6 0 1 8
34 1 0 1 2
34 2 0 1 3
34 3 0 1 9
34 4 0 1 6
34 5 0 1 3
34 6 0 1 5
34 7 1 1 6
35 1 0 1 6
35 2 0 1 6
35 3 0 1 7
35 4 0 1 8
35 5 0 1 9
35 6 0 1 1
35 7 1 1 5
end

Now, I want to fit a discrete-time survival model with the interaction term between the quadratic time variable and a continuous variable.
Could you please show me how to graph the interaction term in Stata code? (M-1sd, M, M+1sd).
Thank you!

logit injury c.time c.safescore c.time#c.time c.time#c.safescore c.time#c.time#c.safescore,or

By the way, what I want is to display the estimated hazard curve showing the effect of safescore while holding all other variables at their mean.

[smith Jason]

When I used margins command, the error message said "factor time not found in list of covariates
r(322);

. margins time, at(safescore = (`mmsd' `m' `mpsd'))
factor time not found in list of covariate

[smith Jason]

As to not using the overall average, the individual values are not independent within ID, and there are an unequal number of observations between IDs. So, some IDs would contribute their presumably correlated values of safescore more than others, making the pooled-value mean representative of Idontknow. In order to help assure that independent observations contribute to the mean, I chose each ID's last observation inasmuch as it is the closest in time to the injury or its most recently observed avoidance.

If you want to get fancy about using some kind of grand mean, then I suppose that you can do something like the following in recognition that the observations are not independent.

mixed safescore i.tim || id: , reml dfmethod(kroger) nolrtest nolog

tempname mean SD
estat sd
scalar define `SD' = r(b)[1, "Residual:sd(e)"]

margins , atmeans df(`e(ddf_m)')
scalar define `mean' = r(table)["b", "_cons"]

local mmsd = `mean' - `SD'
local mpsd = `mean' + `SD'
local m = `mean'

You'd use those local macro variables in place of the ones from yesterday.

But, the purpose is to display the relationship between variations in safescore and variations in hazard, and for that you can just as easily choose three values of safescore arbitrarily. The values of safescore in your dataset range only from six to nine; you could have just chosen the minimum, maximum and a midway value to illustrate the point you want to make on the hazard plot.

As to centering, I'm not sure of what your purpose is in wanting to center the values. Or what would you use to center them (grand mean of whatever flavor?, each individual ID’s mean?).

Can you explain the code below you wrote?
tempname mean SD estat sd scalar define `SD' = r(b)[1, "Residual:sd(e)"] margins , atmeans df(`e(ddf_m)') scalar define `mean' = r(table)["b", "_cons"] Thank you!

[smith Jason]

Can someone help?