Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I have $389$ data points defined in an EventData object and want to judge the quality of fit of this data to a fitted $3$-parameter Weibull distribution. This particular data set contains no censoring but in general there will be censoring in my future data sets. My approach is to find the three parameters using EstimatedDistribution and then to represent the found distribution by generating a number of random variates from it so that they can be compared to the original EventData object. I’m using LogRankTest to produce the $p$-value to judge the fit. My fundamental question is, is there a more elegant way of doing this in Mathematica? My data and code are below.

Here are the raw data:

data1=Table[40.8,{10}];
data2=Table[42.075,{23}];
data3=Table[43.35,{48}];
data4=Table[44.625,{80}];
data5=Table[45.9,{63}];
data6=Table[47.175,{65}];
data7=Table[48.45,{47}];
data8=Table[49.725,{33}];
data9=Table[51,{14}];
data10=Table[53.55,{6}];
data=Flatten[{data1,data2,data3,data4,data5,data6,data7,data8,data9,data10}];

Here is where I define the EventData object and determine the fitted distribution, distribution.

censoring=Table[0,{Length[data]}];
eventdata=EventData[data,censoring]
dist=EstimatedDistribution[eventdata,WeibullDistribution[a,b,c]]

EventData[ ]

WeibullDistribution[2.62045,7.14193,39.7656]

And here is where I non-parametrically determine the $p$-values via LogRankTest.

n=200000;
i=0;
Do[
 {i=i+1,
 If[i>4,Break[]],
 samplevalues=RandomVariate[dist,{n}],
 logrankpvalue=LogRankTest[{eventdata,samplevalues}] ,
 Print["i=",i,"    Log Rank p Value=",logrankpvalue]},
 {5}]

i=1    Log Rank p Value=0.999826

i=2    Log Rank p Value=0.973828

i=3    Log Rank p Value=0.948566

i=4    Log Rank p Value=0.976803

Note that with $n$ set at a relatively large value ($200\,000$), $p$-values of around $0.97$ are produced indicating that the hazard rates of the two data sets are in pretty good agreement. I think this tells me that EstimatedDistribution did a good job in the fitting. Conversely, in general, small values of $n$ (e.g., $10$) produce a much poorer agreement as expected and as evidenced by smaller $p$-values. I generate four $p$-values because I’ve observed that depending upon the particular random values chosen that it is possible to produce a large $p$-value with small values of $n$.

My questions are

  1. Is there a more elegant way of determining the degree of similarity between the EventData object and the found distribution? Perhaps comparing to the theoretical distribution, dist directly instead of sampling from dist?

  2. If this approach (or any improved suggested approach) is robust for the no-censoring case is it also robust if the EventData object contains censored values?

share|improve this question
add comment

1 Answer

One function you might want to look at is DistributionFitTest. Applying this to your data and the estimated Weibull distribution:

DistributionFitTest[data, 
    WeibullDistribution[2.62045, 7.14193, 39.7656],{"TestDataTable", All}]

enter image description here

The values in this table show three different tests for the quality of the fit of the distribution to the data. I'm no expert in statistics, but the fact that these distribution tests seem to show one thing and the LogRankTest seems to show something else might make me nervous. Look closely at what the LogRankTest is doing:

enter image description here

This is not really the same thing as determining the degree of similarity between the data and the estimated distribution, as the OP suggests.

share|improve this answer
    
Bill, Thank you for your inputs. Yes, I first considered DistributionFitTest but ruled it out because I don't think it can handle EventData objects. Moreover, when we look at the p values that you show using DistributionFitTest (and that I confirmed) they are very small ranging from 1E-2 to 3E-158. According to the help documentation for DistributionFitTest "A small p-value suggests that it is unlikely that data came from dist." Therefore how do we reconcile that DistributionFitTest is telling us we have a POOR fit while LogRankTest tells us we have a GOOD fit with p values close to one? –  Steve Oct 12 '13 at 17:03
    
Bill, I agree. LogRankTest and DistributionFitTest are doing two different (subtle to me) things. Neither am I a Statistician but hopefully one will comment. Another reason I favor LogRankTest is that my data is reliability data and I'm actually interested in whether the hazard rates are equal. And yet another reason is that I'm interested in how the populations compare as opposed to only samples of the populations. So I've got lots of reasons to favor LogRankTest. I'd just like to know if there is a more statistically elegant way of doing this than my (crude) "Do Loop Sampling Method". –  Steve Oct 12 '13 at 19:05
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.