Phil Lucht Math & Physics Archive
Home / Math and Physics Files / Math / Probability and Statistics / Survival Analysis

expo model stroke

DOCX · 1.2 MB
Open DOCX file

Personal working notes dated 2/8/15 and 2/9/15, written while Phil worked through the exponential survival model in Lee's survival analysis book. They cover S(t), F(t) and f(t), straight-line fits, right censoring, maximum likelihood for lambda, and chi-square confidence intervals. He applies them to why the CHADS2 study (Gage) error bars differ from his own, noting the notes are rough.

AI-written summary; may contain errors.

Extracted text (machine-read; may contain errors)
Survival Model PhL 2.8.15 Motivation: I noticed that my CHADS2 error bars (from weighted coin flip treatment) were larger than those published by Gage, and he had a footnote that they used an expo survival model, so I studied up a bit on that. I finally figured that they used the second simple model in the Lee book (right here) but still the error bars came out too large. It has to do with how they cut off the data research for each patient, the t2 time. It is an interesting subject associated with the term Poisson Process. These notes are very raw, my first time ever looking at this stuff. I really should write this up better. I also have some PDF's on this subject, Lee is an entire book. First, here is the opening stuff from Lee, Suppose the rate of stroke is determined by something like an exponential survival model, where there is some f(t) which determines chance of having a stroke over time. And example would be S(t) = e-λt You then do a study where you have N patients and you run for T years and you have some subset of the patients n who have strokes, and for each one you have a time ti when they had their stroke. How would you determine λ from such an experiment? Well, the wiki page http://en.wikipedia.org/wiki/Exponential_distribution has some info on this subject. The above model is a constant failure rate, and is memoryless. Here is a good fact: In the very last, we are just applying our general theory of CI to the parameter λ. What do the various Lee functions mean? The rate of having a stroke is a constant over time, that is what λ is. We have S(t) = e-λt dS(t) = -λ S(t)dt dS(t)/S(t) = - λ dt = constant over time This says that in any given time interval dt, the same fraction of people die and this fraction is λdt. As time moves on, the actual number of people who die gets less, because there are fewer people still around and the fraction is constant. Here are the three functions of interest defined above: (all are normalized to a population of 1_ f(t) = λe-λt absolute rate of dying at time t f(t) = dF/dt F(t) = 1 - e-λt cumulative dead at time t S(t) = e-λt cumulative survival at time t. Over time period t, here is what happens probability of patient getting stroke during interval t = S(t)-S(0) = e-λt - 1. ≡ F(t) If λ is very small, this probability is very small. OK, if you start with N patients, then n = N [ e-λt - 1] = expected number of strokes rate over time t = n/N = e-λt - 1 r(t) = e-λt - 1 = F(t) Now maybe I have something! I am at least relating the stroke rate over period T to λ. λ is the instantaneous fractional rate of stroke (so many percent per day) r is the rate of stroke over time interval t Now something that would be easy to measure for a group of N patients would be the cumulative survival situation. Consider Here we start with N patients, and each time 1 gets a stroke, our set of surviving patients downticks by 1. The above red graph then captures the information of WHEN that patient had his stroke. Our task would then be to try to fit the red curve with an exponential blue curve of some λ, and in that way we can compute a value for λ. We are moving forward I think! We could put little dots at the value of N(t) just before each stroke to get this Then we could use just those dots to try and fit the exponential blue curve. I see where logs would help. We have NS(t) = N e-λt ln NS(t) = ln N - λ t or more easily S(t) = e-λt log S(t) = -λt log e Then we could make a linear version of the above plot. Note added later: I found in a pdf (saved) this graph supporting my interpretation of things above, The jaggy line in yellow is the Kaplan-Meier thing, they call KM. For this example, the expo model is not bad, shown in green. The expo is really expo on this linear vertical axis plot, but looks straight. Now first of all, just the fact that the dots all lie on a line supports the expo model. Secondly, we can do a simple fit to obtain λ ! Very simple indeed. I think that is what is being plotted in Example 6.1 from the Lee book. All 25 mice in each test eventually die, but he only has to plot the deaths of the first few to establish that things are linear. The lines don't meet at the origin because this example is also using a second parameter G which is the "guarantee time" where mouse deaths are ignored. Each set of mice is fitted to its own G and λ parameters. So far so good. Now that can be said about this "fitting to the line" situation? This is an Ng like situation. We just try to minimize the variance as a function of λ. (Ng called this variance J). Now we could do this all a different way as follows using the cumulative non-survivor function F(t) We have F(t) = e-λt - 1 so then log F(t) = -(λ loge) t and this would result in a similar straight line test and then fit problem. Let's now jump to Lee Section 7.2 which seems relevant page 176 starts things off, but I am missing what right-censored means so push stack to find that out. OK, section 1.2 gives a very clear description. Censored data means that you don't have all the data! If some patients are alive after the trial stops, you don't know when they would have died. They "go off the end" of the study. In our studies we always know when the event occurred, so there is no "left censoring" there is only "right censoring". The other reason for right censoring is a patient does of something else or withdraws from the study. These are all practical things that really happen. So pop stack and back to Section Section 7.1.1 The idea is to put a + sign after a censored time, fine. Now we have things depending on an Ng vector of parameters! This vector here is called b. Unfortunately, Lee is now going to talk about likelihood functions, and this is a topic that Ng did not address and so I have nothing in my notes. Too bad! He writes where this is a joint probability thing and the first factor group describe people that died at times ti and the second factor group describes things for those who went off the right end and survived. Interesting that once again I encounter the issue for which I am reading the Buck book! So I get the idea here of maximizing this function, but why is it the right function to maximize? Ng was only interested in minimizing his cost function J(θ) to get the best fit for some vector θ. OK, I think I have it. Suppose your two people die at t1 and t2. The probability or "likelihood" of this happening is proportional to f(t1)f(t2), a joint probability in the case of independence. Since this combination of the two deaths DID happen, at those two times, it is fairly probable that those two times caused the probability of this happening be at a maximum. As you vary the parameter λ, you may find that some value of λ causes this event which DID happen to in fact do so at a probability maximum. So consider g(λ) = f(t1;λ)f(t2;λ) = λe-λt1 * λe-λt2 = λ2 e-λ(t1+t2). = λ2e-aλ Now ∂λg(λ) = λ2(-a)e-aλ + 2λ e-aλ = 0 => -aλ2 + 2λ = 0 -aλ + 2 = 0 λ = 2/a Then the answer is that λ = 2/(t1+t2) which agrees with Lee for the general case of n people dying, This solution for finding λ assumes there was no censoring. This means all n patients in the study are "followed to death". Lee goes on to say and we find that the solution 1/λ is simply the average death time of all the people, very simple. So in the case that you can following everyone to death, if you assume an expo survival function, then this gives a very trivial way to find the solution λ. He then goes on to show that wiki result, Notice that (7.2.7) gives an asymmetric interval since the two χ2 objects are not the same. When n > 25 you get the symmetric interval shown. Note that n is the number of people who died and also equals the number of people in this study (they all died). I suspect that if you have a low death rate, the people who did not die don't really matter (they will all have S ≈1), so the analysis above would then apply to our CHADS case where 23 of 512 died. Lee goes on to do a little example of the above, where 21 subjects all die in the study. Using the method above he gets solution λ = .106 per week. And he finds that Is this dramatically asymmetric? 24.4/42 = 0.58 .58 * .106 = 0.061 59.3/42 = 1.41 1.41 * .106 = 0.149 Well, you are 42% below on the left and you are 41% high on the right, so not that asymmetrical. But when converted to a range you get .106 - .062 = 0.044 this is the low side distance .150-.106 = 0.044 this is the high side distance so it is still quite symmetrical, nothing like what I see in the CHADS study. But maybe it gets asymmetric if you convert from λ to rate. Now in fact Lee goes on to do my case where some of the patients survive off the end of the study, and only r of the n people die, so n-r survive. In his example there are 10 people and 5 die and 5 survive off the end. He finds So the difference is that for all those that did not die, you have an extra term in the denom which in his example case is just (n-r) t+. Numerically he then gets So in THIS case we have our 58 sandwiched in the range (19,119) so left width = 58 - 19 = 39 right width = 119-58 = 61 and things are more seriously asymmetric. But we suspect the reason is that r << 25. It is now 9:30 PM on 2/8/15. I have learned a lot about the "exponential survival model", and have seen how you might bet an asym range for λ when the number of deaths is < 25 or so. But the sym ranges persist for larger counts. Recall His largest number of deaths is 25, and for that he has 4.6, 5.9, 7.3 left = 5.9-4.6 = 1.3 right = 7.3 - 5.9 = 1.4 so fairly symmetric. In general the CI ranges are only slightly asymmetric. My puzzle on starting this was merely the question of why they were not perfectly symmetric. Now I know the answer to that question. Feb 9, 2015 First question: Now is the "rate" in the chads paper related to λ? Some definitions: rate1: If you have 500 starting a year only trial and 25 have a stroke, rate1 = 25/500 = 5.00% That is really the rate I have historically been dealing with. Details about the CHADS2 study How do you suppose the "trial" was "run" for chads? It was based on Medicare claims data, so retrospective I guess is the word. Yes, just checked that. A historical study = retrospective study. So how did they do this? Patients are supposed to be AFIB patients. When did each patient "start" in this simulated study? And when did they "end"? Data was gathered in 7 states by certain qualitative people. The first task was to ID those patients who had AFIB, fine. This was determined at the "index hospitalization". At that index event, the properties of each patient were determined, such as high blood pressure. Methods were normalized so all states in the study were done the same way. MEDPAR is a Medicare database which tracks stays at hospitals and nursing facilities. They got dates of death from other sources. OK, but not yet clear to me how this data is used. I understand that for a patient, there is some t1 of an index visit which determines the patient is AFIB and gathers the parameters. But how is t2 determined, the time of a stroke if there is a stroke? Well, first, after excluding various people, they ended up with 1733 AFIB patients each I guess with a start time at the index visit, all in history of course. OK, the t2 time is hospitalization for ischemic stroke. Patients were "censored" either if they died of something else, or if they lived beyond 1000 days of the index visit! Only an initial stroke counted for those with several. So OK, there is how they got t2 times. So here is what we have: (the red bar shows 1000 days) Presumably the study did not include any patients whose index t1 time was less than 1000 days from the time of writing this paper. They probably ruled those out. Then every patient in the study ends either in stroke, other death, or 1000 days duration. That is an interesting way to do it. They then studied the λ parameter, which in this model is called "the hazard rate". It is the instantaneous fractional rate of stroke. They verified the straight line thing! They uses SAS programs called LIFEREG and LIFETEST. They claim to have obtained CI from the "binomial approximation". Stroke rate was based on "patient years of follow-up data". Now it says the 1733 patients were "followed up" for a mean of 1.2 years (median 1.0). What exactly does that mean? Is that after time t2 or in the gap between t1 and t2 ? I think it is time in the gap. But if most patients went 1000 days = 3 years with no stroke, why would the mean follow up only be 1 year? Well it says that during this time there was "readmission for ischemic stroke", so perhaps that means the follow up period starts at time t2 of the original ischemic stroke. This would only serve to determine how many had multiple strokes, or how many died in follow up. OK, I will assume the follow-up starts at t2 just so they can see what happened to the person who had an ischemic stroke. Nothing is said about hemorrhagic strokes! note: 2121/1733 = 1,24 years / patient Now they say there were 2121 patient years of follow up during which 94 people were "readmitted" for ischemic stroke. Of these 94, 71 had one set of codes, and the other 24 had TIA type codes, and they counted these all the same. Now consider: 2+17+23+25+19+6+2 = 94, so the "no of strokes" column in their table does add up to 94. So again I am confused by when the follow up period starts. 27% of the 94 stroke people died in 30 days! They excluded people out of the (65,95) range! They excluded those taking warfarin! This discharge would be from the index hospital visit. The list off the codes of interest. They had a minimum of 1 year of "follow-up" , so I guess yes, follow-up means period after the t1 index visit. So the follow-up range is (365,1000) days, that is very clear. This means follow-up claims in Medicare. surely a fat book like Lee. What does this say? Suppose someone had 2 strokes. They then excluded the 2nd stroke (event), and they excluded patient-days of follow-up after the first stroke. So this sounds like these people were just censured out of the whole study. Now back to that idea that the mean follow-up was 1.2 years. Are they saying that they artificially limited follow-up at random over the group so this was the mean, just to reduce their work required? Maybe this was enough to elicit the results they wanted. Remember that most patients had no stroke at all. Maybe the records did not exist for many of the no-stroke patients. If they never had a stroke admission, there would be no record of anything. So I guess that would not count as follow-up if there was nothing to follow-up. Suppose I align all the t1 times to t = 0 to get this situation: The bottom bar shows 1000 days. The bottom 2 black patients went 1000 days with no stroke, and most of the patients were like this. The top three patients had strokes at the end of the black line segment. They are claiming that the average length of these 3 bars was 1.2 years ? I would expect the average length of these bars to be 500 days since stroke really could happen at any time, and that is 500/365 = 1.37 years, so not too far from 1.2 years. Maybe lots of patients die of other reasons in the first 500 days and are not there for the second 500 days. I suspect that is it. Now, if you add up the lengths of the top 3 bars, they got 2121 patient-years of follow-up, and during that time there were 94 strokes. So this means there are 94 bars that don't reach 1000 days. But this cannot be right because then sum of bar lengths = 2121 years, and only 94 bars, so then each bar would be 2121/94 = 22.56 years long, and we know that is ridiculous. Therefore, "follow-up" must include at least the left end of the full-length bars. So maybe follow-up really was limited in duration by lack of data for all the people. Maybe some people had no records to look at, but that would be in effect follow-up which reveals no stroke. So I now have to assume that "follow-up" includes something more than none of the full length bars, and something less than all of the full length bars. This is very confusing indeed. I need more information about this study to understand what they did. The total length of the full bars is (1733-94)*1000 days which is then (1733-94)*1000/365 = 4490 patient years if none were censured. The total length of the shorter bars had to be less than 94*1000/365 = 257.53 years. So in order to have 2121 years of follow-up, most of that had to be on the full-bars people. What would they mean by "crude stroke rate per 100 patient years" ? 94/2121 = 0.0443 so in some sense for the entire group you get 4.43 % as stroke rate. If you knew how many of the 2121 were in each bin, you might then come up with their crude result. So I guess "crude" means you just divide these numbers. The final adjusted rates are then somehow based on the expo survival model. Adjusted Stroke rate for CHADS2 study. Let's assume they obtained a λ value for each bin of people in the study. How many people then have strokes in a period of 1 year? If nobody got censured during one year, answer would be F(t), which is the total number(as prob) of people who have strokes. So then rate2 = F(1year) where F(t) = 1 - e-λt . So I am claiming rate2 = 1 - e-λ . It would follow that nstrokes = N ( 1 - e-λ), and then annual stroke rate2 = nstrokes/Nstart = ( 1 - e-λ) Now go back to the 94/2121 = 4.43% rate. This says annual stroke rate1 = nstrokes/Npatient years of follow up How are these different? The rate1 seems incomplete since we have this arbitrariness of the amount of follow up, so maybe that is why it is crude. The rate2 comes from the expo model for which you have determined λ,. so rate 2 I think is the thing you really want. So maybe I can indirectly deduce the λ for each bin score 2 bin: annual stroke rate2 = 23/523 = ( 1 - e-λ) Solve to get e-λ = 1 - 23/523 -λ = ln( 1-23/523) = So we find that λ = .04497 = 4.497% If you use the small λ approx you get instead λ ≈ 23/523 = λ ≈ .04398 = 4.398 % Non-stroke deaths were censured so that might affect things a bit. Now let's do this another way using (7.2.10) shown above In my score 2 bin we have r = 23, n = 523, t+ = 1000 days = 1000/365 = 2.74 years (no strokers). The first sum in the denominator has to be less than 23 * 2.74, the second sum is 500 * 2.74, and is likely to be about half that. So I would say λ = 23 / [ 23/2 * 2.74 + 500 * 2.74 ] = So using this method I get λ = 1.64 %/year which is a totally different result and this result is not very affected by the exact size of the first sum. The problem here is that they did not really include all those 1000 days for all the healthy people, so this method fails. If their look window was only 1.2 years on average, maybe a better answer would be larger by ratio 2.74/1.2 and that gives 3.747 %. Now what about CI? In the Lee model which gives where only r people get strokes out of your n people, the error bars are given by : Now in Lee notation, for 95% CI you have α = .05 So for α = .05, Z.025 is the 2.5 % point where 2.5% is in the tail on each side. Now the standard B-1 distribution has σ = 1 and μ = 0 (confirmed on web). So Zα/2 = 1.96 Just by the way, look at his table For 95%, one side must have .4750 of the distribution. This table shows this occurs at z = 1.96. So OK, here is the error formula where r = number two died, and Z = 1.96. The condition for using the above is that n is large, not that r is large! Apply this to our chads = 2 group with r = 23 d = 1.96/ = .4179 Maple 1-d = .5821 maple 1+d = 1.4179 maple Then our range is this left = .5821 * λ left gap = d*λ right gap = d*λ symmetric right = 1.4179 * λ Now in my evaluation above I found for the score 2 bin, λ ≈ .04398 = 4.398 % Then the range is Now recall in our small λ approx we had, annual stroke rate2 = ( 1 - e-λ) ≈ λ = 4.398% If we then simplify identify rate2 with λ, we would claim rate2 = λ = 4.398% CI range = (2.560, 6.236) my computed The data given by the chads study is this for score 2 So they show rate = 4.0 CI range = (3.1, 5.1) they show I am not very close, too bad. My computed range is larger on both sides. Let's try again a different way. We had annual stroke rate2 = nstrokes/Nstart = ( 1 - e-λ) ≈ λ Try λ ≈ rate2f = 4.0 % from the table. Then if r is really 23 we could get which says rate = 4.0 CI range = (2.3, 5.7) my computed Again my computed range is larger on both sides. So I was really unable to use their "hint" that an expo survival model was involved in order to compute their error bars. One more time with full accuracy rate2 = ( 1 - e-λ) = 4.000% e-λ = 1 - .04 = .96 λ = -ln(.96) = .04082 = 4.082 % This says rate = 4.0 CI range = (2.31, 5.52) my computed so the right side is a little smaller. But still my range is much larger than theirs. Conclusion: Due to the total lack of clarity in the chads description of how it modeled this study, I cannot really compute anything at all. They fed all their data into some SAS computer program and out came all these results. Perhaps this kind of thing has been done so many times that professional readers would not question their methods. Here from my recently saved PDF: The list is death times for a set of patients, in monotonic order as Lee does it. Here is the output of the SAS program: Then there is another page So that is rather interesting. Perhaps LIFETEST assumes everyone in the study died and perhaps it uses a simple expo survival model. I think I can get unlimited info on the programs and the theory. What is missing is how the chads people really used their data and this theory. The notion of follow-up is totally hazy. What does chads2vasc have to say on their methods? They also talk about CI using the binomial approximation. They use the ROC curve business. This was a review of past studies, not a hospital records deal. They do give CI ranges, I was wrong about that. Let's skip to has-bled and see what they did. They had 5333 people, study 2003 and 2004. Study described somewhere else. All AFIB people. Did 1-year follow up to assess bleeding. Did use medical records. Their computer program set is SPSS, not SAS. Sounds like they had for each person an exact 1 year long follow-up. I could look up the Fisher test I suppose. They just report out their results. Why are there no CI estimates? Cannot find a web reason. Ran into Sparctool which has some interesting stuff, gets into all the NOACs. Change subject: there is a connection between the expo model and the Poisson world: Change subject: Can I fit the expo survival model into the universe of p(x) pdf examples like dice, height, stroke? There is a pdf which is f(t) = λe-λt and it has a shape For p(x) = f(t), we have μ = 1/λ and σ = 1/λ (wiki). p(x) is probability of having a stroke at time x. Unlike Bernoulli, which is a coin-toss type p(x), this one is continuous. To get the height pdf, you grab some people and measure their heights. The experiment is a person, the outcome is height. For f(g), you grab some people and see who has strokes when. For rolling one die, there is what I call an underlying p(x) which is p(x) = 1/6 for each face if balanced, and p(x) = something else if weighted. For measuring one height, there is some "underlying" pdf p(x) which is probably some bell curve and which you try to figure out by measuring many heights. For measuring stroke rate with the survival model, f(t) is the underlying p(x). But how to you measure it? If everyone had a stroke sooner or later, you could get 1,000,000 people and you just watch them and write down the time when they have their first stroke. You make a histogram with a horizontal axis of 100 years, say, and bin the strokes and you should get the f(t) above. But if we focus on a time interval during which the rate is roughly constant, such as from 2 to 4 p.m. during work days, the exponential distribution can be used as a good approximate model for the time until the next phone call arrives. Similar caveats apply to the following examples which yield approximately exponentially distributed variables: The time until a radioactive particle decays, or the time between clicks of a geiger counter The time it takes before your next telephone call Wiki has a nice page on the Poisson Distribution. The case above where fλ(t) = λe-λt, but the wiki page talks instead about fλ(k) = λk e-λ/ k! which makes no sense to me at all. Elsewhere Well P(x;μ) ≡ e-μ μx / x! Why does this have mean μ ? Why is the variance μ? What is the range of x? Well, x is supposed to be an integer, and it ranges from 0,1,2...∞. So we expect that 1 = ?? = e-μ Σn=0∞ μn/n! = e-μ eμ = 1, fine. μ = E(n) = Σn=0∞ Pμ(n) [n] = Σn=1∞ Pμ(n) [n] = e-μ Σn=1∞ μn/n! * n = μ e-μ Σn=1∞ μn-1/(n-1)! = μ e-μ Σm=0∞ μm/(m)! = μ e-μeμ = μ. E(n2) = Σn=0∞ Pμ(n) [n2] = e-μ Σn=1∞ μn/n! * n2 = μ e-μ Σm=0∞ μm/(m)! * (m+1) = μ e-μ [ μeμ + eμ ] = μ [ μ + 1] = μ2 + μ var = μ2+μ - μ2 = μ there it is. Wolfram shows this Poisson is the limit of the binomial where you set μ = Np and you keep μ a constant as you let N→ ∞. The binomial mean is Np, so it has μ as a mean as well. The binomial is the probability of n heads versus n in a weighted coin toss repeated N times So what is the connection between this Poisson distribution and the expo survival model? I am not sure there is any connection. I think it is the Poisson process that is of interest here. Well there is a connection, but I think I will let it lie. Now go back to fλ(t) = λe-λt as our underlying pdf which has μ = 1/λ and σ = 1/λ. But Lee showed in the simple case that μ = which is just the average of all the times from your experiment, similar to the height study where μ = . So what does our little formula say, σmean = σ / σ = 1/λ and the mean is 1/λ, so σmean = 1/λ / = / Then range of mean would be mean left = - 1.96 / mean right = + 1.96 / and write this as mean left = ( 1 - 1.96 /) mean = mean right = ( 1 + 1.96 /) and this agrees with the earlier quote for this simple case where you invert the thing. So if you did this study, you outcome would be , the mean time for having a stroke (assuming everyone has a stroke). Death would be better, and then = mean lifetime of a person. So finally I have applied the "general idea" that σmean = σ / to this particular pdf. The problem again with chads2 is that we cannot just use the simple expo model because not everybody eventually gets a stroke. We could use the modified model Lee discusses where n-r people go off the right end never getting a stroke. Lets try again on the Lee second model in the chads2 world. We have at score 2 a set of 523 people and of them 23 have strokes. Our mean time to stroke is then = the following expression: = The mean "survival time" here means the mean "time from t = 0 to get a stroke". Now suppose they really did run everybody off the end at 1000 days = 2.74 years. That is just the length of the study as if they all started at once. Then in the numerator, the second sum is 500*2.74 = 1370 years. The first sum we could estimate as 23* 2.74/2 = 31.5 years (assume strokes evenly spread out since rare events). r = 23, so we get = [ 1370 + 31.5] / 23 = 60.93 years I would then argue that this says each year you will have 1/60.93 patients die which is .0164 or 1.64%. But the number chads2 gives is 4%, not 1.6%. The reason this is off is that this is not the method used by chads2. Remember that the mean "follow up time" was only 1.2 years, and this totally clouds the picture. If the study only ran 1.2 years with the same facts as above, you would get = [600 + 13.8]/23 = 26.68 years, and then rate = 3.7%, closer to their number 4. I have a feeling that their study had tricky censuring of the patients in a highly variable manner. Why didn't they just "follow up" all 1733 patients for 2.74 years, using medical records? I can see that for 94 patients who had strokes, they would stop the follow-up at the stroke time. But that is only 1/20th of the patients. Is it possible that a huge fraction of the CHADS2 NRF cohort terminated in a censured manner? If nobody was censured other than the 2.74 years, the study would have 1733*2.74 = 4748.42 patient-years of follow-up. But in practice they only had 2121 such patient-years of follow-up. Only 94 people had a stroke, so they could have removed only 2.74*94 = 258 years at most from the total of patient years. Remember that the NRAF group had average age of 81 !! So perhaps lots of them died? The Gage paper does not say how many of the 1733 patients died less than 1000 days from their index event. If all the people were 81 years old. If you are 80, tables say you get 9 years life expectancy. So how many would die in 2.7 years? Maybe 2.7/9 = 0.3 = 30%, so that is quite a few. That would knock out 4748* .3 = 1424 years / 2 since they would die half way on average, so you only lose 712 patient years that way, and you end up with about 4000 py but they only had 2121 py that got used. Where is the other half?? Idea: maybe you have to have a positive tracking active to really count a patient year. If you simply don't hear anything on the IDC channel, that does not mean they did not have a stroke. You might have just lost contact and you don't really know anything about that patient, so that patient is censured when it leaves the IDC radar. But how could a Medicare history stop other than death? Idea: maybe there were limited resources for doing the IDC research, and they could not trace every single person all the way over 1000 days. On average, they only tracked for 1.2 years, and that is around half the 2.74 and maybe that is where the missing half went. One idea of the censuring deal is that people "contribute" to the study even if their time in the study is right censured, and the reason for such a right censure might be that it cost too much in $ or time to track all the years. Just because a person is in this way censured by the researchers, they still contribute to the study for the time they were active in it. Maybe different states had different resources, and some could track all the way out to 1000 days, but most could not. So in my group of 523 of whom 23 had strokes, it may be that half the data is lost. So my assumption that 500 of this group made it through 1000 days with no stroke is just wrong. Maybe only 50 of them did that for all I know. Then basically for every patient of the 1733, there was a t1 at their index event, and there was some t2 that was either stroke, death, withdrawal, or "termination" by the research team. Then you have 1733 pairs of the form (t1, t2) where all you really know is that t2 - t1 < 1000 days. Then you have to use the Lee formula, where the t(i) are times of stroke, while t(i)+ are times of censure. I have no access to this information, so I cannot compute the exact from this formula. But maybe I know this much Σi=1r ti + Σi=r+1n ti+ = 2121 Then we know that = 2121/94 = 22.56 years = 1/λ This would then give an average stroke rate for the entire group, rate = 1/22.56 = 0.0443 = 4.43% = λ which is in the right ballpark I would guess. I did an estimate of this perhaps. In my 4.7 I get <rate> = 4.83 which is not too far from the 4.43. Comment: λ is the fraction getting stroke at any time. At t = 0 all the patients are there. If we find that λ = 4.43%, then that is the fraction that are stroking just after the t = 0 start of the race. And that I think is the number of interest. Now go back to earlier : d = 1.96/ = .4179 Maple 1-d = .5821 maple 1+d = 1.4179 maple This seems correct for the score = 2 bin. But the λ value I now think is 4.0 and is the same as the rate they quote. So then rleft = 4.0*.5821 = 2.328 rright = 4.0 * 1.4179 = 5.67 and then this says rate = 4.0 % CI = (2.3%, 5.7%) my compute whereas they give rate = 4.0 CI = (3.1, 5.1) they show so as before, my error bars are larger than theirs on both sides, but OK ballpark. Now I figured out the χ2 table in B-2 and redid the above using that method and I found Lfactor = .63 Rfactor = 1.43 and then rleft = 4.0*.63 = 2.52 rright = 4.0 * 1.443 = 5.772 rate = 4.0 CI = (2.5, 5.8) using χ table so I still don't get much closer to their claim. But I did this wrong, because the Now it is true that the † thing on chads says "assuming aspirin was not taken". How would that affect things? So quite a few were taking aspirin. Maybe the 523 and 23 include aspirin takers! aspiring takers = 529/1733 = 0.31 = 31% Roughly, suppose then we use 69% to knock things down so there were then only 23*.69 = 16 strokes in this group who took no aspirin. But having a smaller r just makes the error bars larger. whereas I would like to have them be smaller. As r gets smaller, the χ22r get smaller, but the denominator r also gets smaller, so not clear what will happen. Lets try r = 16 and then we have χ232,.975 = 18 χ232,.25 = 48 2r = 32 left factor = 18/32 = 0.56 right factor = 48/32 = 1.5 .56*4 = 2.2 1.5*4 = 6.0 rate = 4.0 CI = (2.2, 6.0) using χ table As I expected, smaller r makes the error band larger, just as it does with the large n simple formula. , But in reality, they know exactly how many strokes arose in the aspirin group and I don't know that. But it must be at most 23 so no fiddling with that is going to make smaller error bands. Conclusions: (1) the Lee model is not really applicable to what chads2 is doing (2) the chads2 people have too-small wrong CI numbers (3) I am misunderstanding something