Dear Stata users,
I have been running simulation studies investigating interrupted time series studies.
My original computer used Stata 15.1 and I updated to Stata 16.1.
For reproducibility I have ensured to set versions and save random seeds.
When using mixed with restricted maximum likelihood I noticed that some of the autocorrelation estimates were exactly equal to zero (while the true autocorrelation was 0.8) when using one computer and not when using another.
I later tested this on a computer that had both Stata 15.1 and Stata 16.1 installed and fully updated and found the same issue.
I went back and checked some of the datasets (one of which is at the end of this post, but there were many).
I have made several interesting/confusing/potentially annoying observations.
1. the same dataset with the same do file gives a different result when using Stata 15.1 or Stata 16.1 when both are fully updated, using the same computer or different computers even when version 15.1 is included in the do file.
2. the restricted maximum likelihood default maximizing technique, Newton-Raphson, sometime does not work as expected (giving no estimate of autocorrelation and incorrectly converging), but no indication of this (such as an error message) is shown.
When using the same dataset and do file:
Stata 15.1 gives autocorrelation estimate of .7849656 for both the NR and DFP maximizers.
Stata 16.1 gives autocorrelation estimate of 0 for NR (with no indication of any difficulties) and .7849656 for DFP.
And for further confusion:
When copying the below code into a do file (with the data input from using dataex)and running it:
Both Stata 15.1 and Stata 16.1 give autocorrelation estimates of 0 for NR (with no indication of any difficulties) and .7849656 for DFP.
This doesn't make sense to me, as the data should be the same as the original dataset (which gave 0 in Stata 16.1 for the NR and .7849656 for Stata 15.1)!
I've tried this with several different datasets now and sometimes get different results when copying from the do file into the code window below, then copying from the code window back into Stata. I don't know how/why this could be!?
If anyone could shed some light on any of these issues (getting different results using different Stata versions; REML with NR not giving expected feedback; getting different results when copying/pasting) that would be marvellous.
Please find below:
1. screenshot of the two Stata windows on my PC showing different results from the same do-file and dataset.
2. the code used to produce the data and graphs.
Cheers,
Simon.
I have been running simulation studies investigating interrupted time series studies.
My original computer used Stata 15.1 and I updated to Stata 16.1.
For reproducibility I have ensured to set versions and save random seeds.
When using mixed with restricted maximum likelihood I noticed that some of the autocorrelation estimates were exactly equal to zero (while the true autocorrelation was 0.8) when using one computer and not when using another.
I later tested this on a computer that had both Stata 15.1 and Stata 16.1 installed and fully updated and found the same issue.
I went back and checked some of the datasets (one of which is at the end of this post, but there were many).
I have made several interesting/confusing/potentially annoying observations.
1. the same dataset with the same do file gives a different result when using Stata 15.1 or Stata 16.1 when both are fully updated, using the same computer or different computers even when version 15.1 is included in the do file.
2. the restricted maximum likelihood default maximizing technique, Newton-Raphson, sometime does not work as expected (giving no estimate of autocorrelation and incorrectly converging), but no indication of this (such as an error message) is shown.
When using the same dataset and do file:
Stata 15.1 gives autocorrelation estimate of .7849656 for both the NR and DFP maximizers.
Stata 16.1 gives autocorrelation estimate of 0 for NR (with no indication of any difficulties) and .7849656 for DFP.
And for further confusion:
When copying the below code into a do file (with the data input from using dataex)and running it:
Both Stata 15.1 and Stata 16.1 give autocorrelation estimates of 0 for NR (with no indication of any difficulties) and .7849656 for DFP.
This doesn't make sense to me, as the data should be the same as the original dataset (which gave 0 in Stata 16.1 for the NR and .7849656 for Stata 15.1)!
I've tried this with several different datasets now and sometimes get different results when copying from the do file into the code window below, then copying from the code window back into Stata. I don't know how/why this could be!?
If anyone could shed some light on any of these issues (getting different results using different Stata versions; REML with NR not giving expected feedback; getting different results when copying/pasting) that would be marvellous.
Please find below:
1. screenshot of the two Stata windows on my PC showing different results from the same do-file and dataset.
2. the code used to produce the data and graphs.
Cheers,
Simon.
Code:
* Example generated by -dataex-. For more info, type help dataex
clear
input float(time observation intervention level_change slope_change slope_pre)
1 3.0149505 0 0 0 1
2 3.031883 0 0 0 2
3 3.1242235 0 0 0 3
4 4.2264733 0 0 0 4
5 4.0782666 0 0 0 5
6 3.482328 0 0 0 6
7 1.4031403 0 0 0 7
8 1.4085115 0 0 0 8
9 1.0365456 0 0 0 9
10 -.29633772 0 0 0 10
11 -1.677048 0 0 0 11
12 -1.623045 0 0 0 12
13 -.4066716 0 0 0 13
14 .6809254 0 0 0 14
15 .5246355 0 0 0 15
16 -.5750102 0 0 0 16
17 -1.0658898 0 0 0 17
18 -1.0694761 0 0 0 18
19 -.27314278 0 0 0 19
20 .4933012 0 0 0 20
21 -.5803595 0 0 0 21
22 -.8530504 0 0 0 22
23 .3393627 0 0 0 23
24 -.3632981 0 0 0 24
25 -.7547657 0 0 0 25
26 1.1019006 0 0 0 26
27 1.501023 0 0 0 27
28 .4218852 0 0 0 28
29 2.516184 0 0 0 29
30 1.0772717 0 0 0 30
31 -1.353526 0 0 0 31
32 -1.4707215 0 0 0 32
33 -.20380923 0 0 0 33
34 .9197118 0 0 0 34
35 1.9651598 0 0 0 35
36 1.2351445 0 0 0 36
37 .863229 0 0 0 37
38 -.03003758 0 0 0 38
39 1.0612937 0 0 0 39
40 -.20074587 0 0 0 40
41 -1.6133785 0 0 0 41
42 .7977557 0 0 0 42
43 1.8189863 0 0 0 43
44 .8851298 0 0 0 44
45 -.02789556 0 0 0 45
46 .21301727 0 0 0 46
47 -.07175163 0 0 0 47
48 -.6914082 0 0 0 48
49 .06886576 0 0 0 49
50 -.517148 0 0 0 50
51 1.79529 1 1 0 51
52 .6271598 1 1 1 52
53 1.0698519 1 1 2 53
54 2.39213 1 1 3 54
55 3.759435 1 1 4 55
56 3.714929 1 1 5 56
57 2.2656763 1 1 6 57
58 -.016826201 1 1 7 58
59 2.3457026 1 1 8 59
60 3.610192 1 1 9 60
61 3.473347 1 1 10 61
62 4.3406916 1 1 11 62
63 4.783299 1 1 12 63
64 5.872327 1 1 13 64
65 5.403982 1 1 14 65
66 5.146616 1 1 15 66
67 5.015276 1 1 16 67
68 4.72711 1 1 17 68
69 4.4966965 1 1 18 69
70 4.7581396 1 1 19 70
71 4.43996 1 1 20 71
72 4.0046325 1 1 21 72
73 4.127524 1 1 22 73
74 5.468728 1 1 23 74
75 4.058884 1 1 24 75
76 1.8522925 1 1 25 76
77 2.1346645 1 1 26 77
78 2.701144 1 1 27 78
79 3.684405 1 1 28 79
80 1.9640895 1 1 29 80
81 1.8081708 1 1 30 81
82 2.0514982 1 1 31 82
83 2.830338 1 1 32 83
84 2.05137 1 1 33 84
85 4.111118 1 1 34 85
86 4.288389 1 1 35 86
87 4.140128 1 1 36 87
88 5.226665 1 1 37 88
89 5.798832 1 1 38 89
90 4.707952 1 1 39 90
91 5.062461 1 1 40 91
92 5.913624 1 1 41 92
93 5.432038 1 1 42 93
94 2.3810606 1 1 43 94
95 1.6668195 1 1 44 95
96 2.634732 1 1 45 96
97 4.1317163 1 1 46 97
98 4.184508 1 1 47 98
99 4.2190104 1 1 48 99
100 5.312669 1 1 49 100
end
local intervention_time = 50
// run standard REML (which uses Newton-Raphson)
mixed observation slope_pre level_change slope_change, res(ar 1, t(time)) var reml
predict mixed_obs
matrix local_results_mixed = r(table)
local num_cols = colsof(local_results_mixed)
local rho_est_mixed = tanh(local_results_mixed[1,`num_cols'])
local display_rho_est_mixed : display %3.2f `rho_est_mixed'
// run REML with the dfp alternate
mixed observation slope_pre level_change slope_change, res(ar 1, t(time)) var reml tech(dfp)
predict mixed_dfp_obs
matrix local_results_mixed_dfp = r(table)
local num_cols = colsof(local_results_mixed_dfp)
local rho_est_mixed_dfp = tanh(local_results_mixed_dfp[1,`num_cols'])
local display_rho_est_mixed_dfp : display %3.2f `rho_est_mixed_dfp'
// graph and display the results
twoway scatter observation time, msym(X) mcol(blue%30) ///
|| lfit mixed_obs time if time < `intervention_time' , lcol(blue) || lfit mixed_obs time if time >= `intervention_time', lcol(blue) ///
|| lfit mixed_dfp_obs time if time < `intervention_time' , lcol(orange) || lfit mixed_dfp_obs time if time >= `intervention_time', lcol(orange) ///
graphregion(col(white)) legend(order(1 "values" 2 "NR" 4 "DFP") rows(1) region(col(white))) ///
ylab(,angle(h)) title("{&rho}{sub:true} = 0.8, {&rho}{sub:NR}= `display_rho_est_mixed', {&rho}{sub:DFP}= `display_rho_est_mixed_dfp'")
Attached Files
