How to assess effect of intervention in one state versus another using annual case fatality rate?

I am a beginner in statistics with just basic knowledge. I have these data: cases, deaths and CFR (Case Fatality Rate-deaths per 100 cases) of a disease for 17 years (1994-2010) from 2 neighbouring states where people can walk across the states freely. This is a population based cohort study.

Data are available from 1994. The treatment protocol was started in 1996. State 1 and state 2 implemented the same treatment but state 2 could not implement the treatment perfectly due to local administrative problems. The death rate fell in 1 state and continued to be high in the other state. Because it would be unethical to subject patients to do a case-control study, I want to analyse the following available data to see if there is a significant difference in the death rates and risk ratios of these 2 states because of poor implementation of the treatment guidelines in state 2 from 1996 to 2010.

St. refers to State

Year   St.1 Cases  St.1 Deaths St.1 CFR St.2 Cases St.2 Deaths  State2 CFR  Risk Ratio
1994    1836        383         20.86     583        121         20.75        0.99
1995    1246        257         20.63    1126        227         20.16        0.98
1996    1450        263         18.14     896        179         19.98        1.10
1997    2953        407         13.78     351         76         21.65        1.57
1998    1161        149         12.83     1061       195         18.55        1.43
1999    2924        434         14.84     1371       275         20.06        1.35
2000    1729        169          9.77     1170       253         21.62        2.21
2001   1888      275         14.57     1005      199          19.80        1.36
2002    919     178          19.37      604      133         22.02         1.14
2003    865        142         16.42      1124      237         21.09         1.28
2004    1543       131          8.49      1030      228         22.14         2.61
2005    2887       336         11.64      6061     1500         24.75         2.13
2006    1484       108         7.28       2320     528          22.76         3.13
2007    1592       75          4.71       3024     645          21.33         4.53
2008    1920       53          2.76       3012     537          17.83         6.46
2009    1477       40          2.71       3073     556          18.09         6.68
2010    1534      26           1.69       3540     494          13.95         8.23


Kindly advise me what is the best way to go.
I have basic knowledge of using SPSS V19 and Comprehensive Meta-analysis V2.

-
@mbq: Thank you for your time in editing this –  DrWho Mar 16 '11 at 13:40
@DrWho Is this intended as a replacement for your previous question at stats.stackexchange.com/q/8187/919 ? –  whuber Mar 16 '11 at 16:17
@whuber: I am amazed at your keen observation. Yes. Now we could recover more data from a crashed hard disk. I thought that increase in the sample size will increase the power of the study and yield better inferences. –  DrWho Mar 17 '11 at 0:44
@DrWho It might be best to merge these two questions, then. That can help garner more focused answers for you. Consider editing your previous one in light of the information in this one. –  whuber Mar 17 '11 at 2:05
@whuber: I agree that this question was answered but the new data (17 years ) may present analysis opportunities unavailable with the smaller data set (6 years). IrishStat 2 hours ag –  IrishStat Mar 17 '11 at 13:39
show 6 more comments

I don't know why I took the time to answer this. Is it because I can or maybe it's because DrWho seems to think it is very important. In either case ....

Though well intentioned

“Time series expert modeler of IBM SPSS Forecast v19 was used. Both exponential smoothening models and ARIMA models were examined.Outliers were detected and prevented from influencing parameter estimates”

may have suffered from an inability to detect level shifts i.e. “Intercept Changes” which are a sequence of pulses with the same value and sign. Note below a reasonable model for STATE1 using all 17 values suggests a Level Shift at 2004 ( period 11). This model [AR(1)] was used to cleanse STATE1 of unspecified background factors that may have been present to cause significant changes in Y given X.

Y(T) = -87.899
+[X1(T)][(+ .158)] M_CASES +[X2(T)][(- 77.4369)] :PULSE 7 I~P00007STATE1 +[X3(T)][(- 65.2775)] :PULSE 15 I~P00015STATE1 +[X4(T)][(- 112.77 )] :LEVEL SHIFT 11 I~L00011STATE1 +[X5(T)][(+ 43.0999)] :PULSE 9 I~P00009STATE1 +[X6(T)][(- 58.4117)] :PULSE 4 I~P00004STATE1 + [(1- .840B** 1)]**-1 [A(T)]

Leading to a cleansed set of values FOR STATE1

1994 383.0000000000
1995 257.0000000000
1996 263.0000000000
1997 465.4116693551
1998 149.0000000000
1999 434.0000000000
2000 246.4369202361
2001 275.0000000000
2002 134.9000542626
2003 142.0000000000
2004 131.0000000000
2005 336.0000000000
2006 108.0000000000
2007 75.0000000000
2008 118.2775018144
2009 40.0000000000
2010 26.0000000000

Notice that a simple line graph between Y and X visually support the change in the relationship between Y and on or about period 11 (2004) such that the Y values are clearly lower than expected for the period 2004-2010 ( 11-20) as compared to period 1994-2003 (1-10). This is a classic case of an outside factor impacting either Y or X ( but not both !) starting at time period 11. Normal statistical commentary would refer to this level shift as a “lurking variable” confounding simple analysis if untreated.

For STATE2

Y(T) = -6.4072
+[X1(T)][(+ .213)] M_CASES +[X2(T)][(- 254.25 )] :PULSE 17 I~P00017STATE2 +[X3(T)][(+ 214.32 )] :PULSE 12 I~P00012STATE2 +[X4(T)][(- 92.6947)] :PULSE 16 I~P00016STATE2 +[X5(T)][(- 98.6907)] :PULSE 15 I~P00015STATE2 +[X6(T)][(+ 39.8298)] :PULSE 13 I~P00013STATE2 + [A(T)]

1994 121.0000000000
1995 227.0000000000
1996 179.0000000000
1997 76.0000000000
1998 195.0000000000
1999 275.0000000000
2000 253.0000000000
2001 199.0000000000
2002 133.0000000000
2003 237.0000000000
2004 228.0000000000
2005 1285.6763010529
2006 488.1702389682
2007 645.0000000000
2008 635.6907401611
2009 648.6947149772
2010 748.2497352909

Note that the Pulses were identified GIVEN the number of cases.

Now proceeding with the CHOW TEST to test the significant difference between the two sets of regression coefficients:

We get the following OLS Model for the combined 34 (cleansed values used )

With Error Sum of Squares = 363716.

For STATE1 we get

Y(T) = -49.293
+[X1(T)][(+ .150)] M_CASES
+ [A(T)]

Sum of Squares 142759.

And for STATE2

Y(T) = -6.4072
+[X1(T)][(+ .213)] M_CASES + [A(T)]

Sum of Squares 1445.64

As before we have the following F test

Numerator = [363716 – (142759+1445)] /2 = 109,776

Denominator = 363716/30 = 12,123

Computed F value of 9.0 with 2 and 30 degrees of freedom is significant at alpha less than .001 . Thus one could conclude that there is a statistically significant difference between the two states at about a 99.9% Level of Confidence given that there was a significant effect in STATE1 at or about 2004.

-
@IrishStat: Thank you for your time and expertise. I can see how much valuable time you have spent to do the calculations. Since I am a beginner in Statistics, I require some time to understand your calculations and method. I will get back at the earliest. Kind Regards –  DrWho Mar 21 '11 at 3:15
@IrishStat: I am extremely happy with your excellent guidance. It is taking a long time for me to understand your calculations. For state 1, the model was ARIMA (1,0,0). For state 2 what was the model? Is it the same? Can you kindly let me know? –  DrWho Mar 21 '11 at 8:20
@IrishStat: While trying to search for more information on OLS model, I chanced upon your very useful document autobox.com/pdfs/Pickett.pdf published in THE JOURNAL OF BUSINESS FORECASTING, SPRING 2005. Now I know why it is taking so much time for me to understand what you wrote in a few minutes. i am so happy my call for help through this website attracted the attention of a great and dedicated man. –  DrWho Mar 21 '11 at 8:48
@DrWho: The model for STATE1 is a Transfer Function Model because it includes causal variables ( # of cases plus the Interventions that were Detected ) and an ARIMA component (1,0,0). The model for STATE2 was similar but had an ARIMA component (0,0,0) which means it was a simple multiple regression model (OLS) . Regression(OLS) is a particular case of a Transfer Function. –  IrishStat Mar 21 '11 at 12:43
@IrishStat: Thank you very much. Every time you reply, I learn something more to my great advantage. I am still struggling to understand the meaning of your approach. Best Regards –  DrWho Mar 21 '11 at 12:47
show 2 more comments

First, you need to decide what questions you want to ask. Are you comparing the states, year by year? Are you looking at changes within state by year? Are you comparing the three variables, within states? Are you doing something else?

To help make this decision I would make a line graph like @CaseyTsui suggested. As Yogi Berra said "you can see a lot by looking". But I would make three separate graphs (one for each variable), each with two lines (one for each state). I would also make two graphs (one for each state) with three lines on each (one for each variable) after rescaling so that they are on similar scales.

If the graphs don't answer your question, then, as @IrishStat said, you need to look at time series analysis. This is a complex field, but a good first step is to look at the autocorrelation function for each state; this can be done in SAS/ETS or in R (and very likely in other packages as well, but those are the two I know). Then you need to start looking at things like ARIMA models; the material in the documentation for PROC ARIMA in SAS is pretty good, and not too complex.

But that may not be necessary at all. Just eyeballing the data tells me that there are big differences. Don't get needlessly formal. These graphs may pass the IOTT; the interocular trauma test - it hits you between the eyes. You could make many other plots, depending on what you are interested in. Then post the plot and say "Behold!"

(Of course that won't give you things like p-values; but (to paraphrase my favorite stats professor, you don't need to p on the research).

-
If you ignore cases & deaths and simply study rates then ARIMA might be relevant. The question is "Are the coefficients for the appropriate ARIMA Model the same across the two groups(states). One answers this by estimating the ARIMA model coefficients GLOBALLY and then SEPARATELY in order to construct the F test.In general one would not want to ignore the two observed series(cases and deaths) but to construct an ARMAX Model and perform theTest on the ARMAX coefficients.Care would have to be taken to deal with Pulses,Level Shifts,Seasonal Pulses/Local time trends in either case. –  IrishStat Jun 22 '11 at 14:49
Though results are impressive without any statistical analysis, applying correct statistical method has been a problem for me. This looks to be a complex issue. I am eagerly waiting for a complete and step-wise advice. –  DrWho Jun 23 '11 at 2:56
@DrWho A) You didn't answer my questions. B) Yes, it is a complex issue. It is therefore probably not wise to expect a complete and stepwise solution posted on the web. This is a question for which you really need to hire a consultant (someone who knows time series) and work with them over a fairly long period to figure out what you want to do. You can get some advice here (I have certainly gotten good advice here!) but in this case, you need more than what will fit. –  Peter Flom Jun 24 '11 at 10:33
Sorry for the delay. Thank you for your query. My aim is to show that the different treatments in 2 states made a difference in the outcome of cases. –  DrWho Jun 25 '11 at 7:37

What time-lag might you expect between cases recorded, and fatality? What time-lag between start of treatment and impact on fatality rates?

If either of those numbers is much greater than one year, then there may be a case for aggregating all your data from first year of treatment impact (i.e. 1996+time to impact) to 2010, and just test to see if mean CFR rates vary. If possible, ask whoever rejected the time-series approach whether this would satisfactorily deal with the confounding factors that concern them.

Do look at specifics of the confounding factors: for example, if palliative care for terminal patients was better in State 2, then those patients originally in State 1 for whom treatment had failed, might move to State 2 for their final days/weeks/months of life; in that case, their deaths would be registered in state 2, and the numbers you have won't provide you much useful information, unless you could get the numbers of terminal migrations between states.

-
Thank you for the stimulating questions. 1. Time lag between case recording to death is up to 6 months. Cases recorded in 1 state are followed up for 6 months. No case was ever lost to follow up. Case recorded in December 2007 and dying in May 2008 was recorded as 2007 case only by the state. 2. When a journal refuses a paper with a comment, they don't give us a chance to explain. So I want to send to another journal. –  DrWho Mar 20 '11 at 9:55
3. The native place of the patient and the detials of prior treatment were recorded at admission it self. This work was done with the help of an International Organizations. Migrated cases were taken as belonging to the state of their origin. 4. I would be extremely happy to clarify any other point. –  DrWho Mar 20 '11 at 9:55
I wrote like this. Selection of the model for forecast generation: Time series expert modeler of IBM SPSS Forecast v19 was used. Both exponential smoothening models and ARIMA models were examined. Outliers were detected and prevented from influencing parameter estimates. Maximum number of Lags in ACF and PACF output checked were 24. Stationery R-square (larger values indicate better fit) was the goodness of fit measure for selecting the best fitting forecasting model. –  DrWho Mar 20 '11 at 10:06
CFR from 1994 to 2005 were used for forecasting CFR along with upper and lower confidence limits from 2006 to 2010 and compared with the actual observed values till 2010 and the results were interpreted. Significance at 99% confidence level and objective endpoint (death) were used to minimize the effect of the confounding variables and make the conclusions more reliable. The data from the same state were compared. Multiple areas-India, Both states were studied independently and conclusions were drawn. Please tell me if my methodology was wrong. –  DrWho Mar 20 '11 at 10:07
Please ignore Multiple areas-India in the above comment.That came in by mistake –  DrWho Mar 20 '11 at 10:18

To represent the data visually, you can do a simple line graph:

x-axis: year

y-axis: CFR

Stratify by state.

For a test statistic determining whether the CFRs for each state are significantly different from each other over time, you could do an ANOVA between Year, CFR, and State as the third variable. You'd have to first reshape the data to long format. The significance of the State F-statistic will let you know if there is a unique difference in the CFRs between the two states.

-
Casey: This is time series data thus the observations are correlated thus you can not use the ANOVA approach as you don't have independent observations –  IrishStat Jun 21 '11 at 22:27
@IrishStat: Of course, you are right. –  Casey Tsui Jun 21 '11 at 23:26
I am a beginner in statistics. I am still waiting for the right method. Kindly guide me step-wise. –  DrWho Jun 22 '11 at 4:06

With only just these two cases, you cannot reliably estimate a treatment effect, but you can summarize your data as follows, assuming that the the number of deaths is a draw from a Poisson distribution.

\begin{align} &\Pr(\text{Deaths}_{it}) = \lambda_{it}(\text{Cases}_{it}) \\ &\ln(\lambda_{it}) = \beta_{0i} + \beta_{1i}t + \beta_{2i}t^2 + \beta_{3i}t^3 + \ldots \\ &\beta_{0i} = \gamma_{00} + \gamma_{01}\text{State}_i \\ &\beta_{1i} = \gamma_{10} + \gamma_{11}\text{State}_i \\ &\beta_{2i} = \gamma_{20} + \gamma_{21}\text{State}_i \\ &\beta_{3i} = \gamma_{30} + \gamma_{31}\text{State}_i \\ &\ldots \end{align}

Where $\text{Deaths}_{it}$ are the number of deaths in $\text{State}_i$ in $\text{year}_t$; $\lambda_{it}$ is a state and year specific rate of death, and $\text{Cases}_{it}$ are the number of cases in a state in a given year; $t$ is the year of observation; and $\text{State}_i$ is an indicator variable that takes on a value of 0 for the first state and 1 for the second state.

Using a Poisson regression in a general linear models package you can estimate the state specific differences in the average decline ($\gamma_{10}$ and $\gamma_{10} + \gamma_{11}$) and "acceleration" (the effect of the higher order $t$ terms), if you use orthogonal polynomials for $t, t^2, t^3, ...$.

In R, you can use poly and the glm functions:

glm(deaths ~ (0 + cases + cases:poly(year, 4)
+ state:cases + state:cases:poly(year, 4)),
family="poisson",
data=cfr))


With this data, including up to a quadratic term seems to amply summarize any trend.

-
Thank you for your interest in guiding me. I will try this and get back –  DrWho Jun 24 '11 at 4:12