# What R packages or algorithms would you use in that kind of analysis

I've tried to make a short summary of what I'm trying to do but it wasn't clear at all, so sorry for the length of the post.

I'm currently trying to study the length of sick leaves in administrations and so i have a database where have been put the lengths, dates, and other infos for thousands of leaves.

I'm a complete beginner in that type of analysis. First of all, I'm studying the percentage of sick days that happens beyond the x-th day of a sick leave*. I've used the survival analysis package (survival) to simplify my work but this package is focused specifically on survival analysis and the ratio that i study isn't the survival ratio. So i'm looking for a package or a method that could help with what i'm actually looking at.

*The actual ratio that I'm studying is the number of sick leave's days that happened after the x-th day divided by the total number of sick leave's days. (For instance if i have 3 leaves that are 2, 8 and 10 days long : the ratio for day 5 is (3+5)/(2+8+10)=0,4)

The data that i have is a dataframe where in each line I have a single sick leave with pretty much any information that could be useful in a statistical analysis BUT I'm not supposed to use every piece of information. I'm supposed to study this in order to help for the pricing of contracts and the only variables that are communicated to them (and therefore are relevant in my work) are :

-the type of contract (there are two types : I(mostly part-time workers) or C(full-time)),

-the type of sick leaves (3 types : Ordinary(short <1 year), Rare(<3 year), very Life threatening(<5years))

-I also have the amount of deductible days in the contract of the person who had a sick day (to know wheter or not deductible days induce a bias)

-A few other things (I'll go into details if it is necessary)

-Also my data was censored so I have a column with 0s for censored data and 1s for data that isn't.

The only thing that has been taken into account for the censoring is this : the last time the dataframe has been updated was the 31/12/2013 any sick leaves that was still going on at that time has been considered censored. To be more accurate, what I have is the date on the medical report that allowed the sick leave so for instance I can have a sick leave where the end date is 25/02/2014 but it is not the actual end of the sick leave cause afterward it might have been extended so I put it has been censored at the 25/02/2014.

First of all in order to correct censoring I used the survival package like so :

  s <- Surv(time2013,event2013)
fKM <- survfit(s ~ 1,data=mydata)
res <- summary( fKM)


Afterwards I only considered the results that were given to me by the survival model.

What I've been asked to do was to determine what variable had enough influence. (I haven't been told what enough actually was so what I did so far, was to say if, for instance, when I compare sick leaves where the work contract was : I(part-time) with C(full-time), if I have more than ~4-5% of difference in my ratios then the variable had an influence and so I divided the population in two categories that i studied separately.)

The problem I have now is to find a statistical way to tell if a variable has enough influence, so I have the following plot (the two types of contract) I've said that there was enough difference so that when pricing the contracts they should price their insurance differently for each type of contract. Red is I(part time), Black is C(full time) the plot tells me that full-time workers have slightly longer ordinary(<1 year) leaves.

What i would like to do now is find a way(most likely with an algorithm or another package) to measure the difference between those two plots. In order, to be able for instance to compare this difference with the difference there is in rare leaves(<3 years). So that I can know if a difference is significant by comparing it to the other differences.

If you think another approach would be more useful I have no problem with this, comparing the differences was just what seemed the most intuitive analysis to me.

Considering the Cox model, unless I'm making an error, I don't think (I might be wrong) that it is working in my situation because of the results that I get.

For instance I have two situations ; the first one here is the plot of the ratio I'm studying (not the ratio that the cox model is using, I didn't want to put too many images the post is already too long) for a population that is divided in 2 groups (Red is full-time workers and Black is part-time workers). I apply "coxph" I get a coefficient for the type of contract of 0.223 which is interpreted, if I understood correctly, as : "full-time workers have longer sick leaves", which we don't really see on the plot of the percentages of days beyond the x-th day :

But on another situation where I have this plot :

When I use "coxph" I get a coefficient of 0.2175 which is basically, giving me the same conlusion. But intuitively i would have said that in the first situation the differences were insignificant and in the second they were meaningful. That's why I haven't used "coxph". Maybe you can find an issue in the code that i used :

mydata_MO$SI_DUREE_ARR is the length of the sick leaves; mydata_MO$NON_CENSURE is 0 for censored data and 1 for not censored data;

IND_AG_CONTRAT is 1 for full-time, 0 for part-time

timeMO<-mydata_MO$$SI_DUREE_ARR eventMO<-mydata_MO$$NON_CENSURE
coxph.fit <- coxph(Surv(timeMO,eventMO) ~ IND_AG_CONTRAT, data=mydata_MO)
summary(coxph.fit)


Sadly, I can not share the data, but what i believe is the issue is that both cases are very similar except that there are a few extremely long leaves in the second case that are not in the first situation which is why there is such a big difference on the graphs.

• While the Cox model (using pkg:survival with coxph) would give you access to all the information you needed, you will need more than just the coefficients. The curves make me think that a more parametric analysis would be justified (also supported by pkg:survival). The knowledge needed to answer it would be that gained in a second year graduate course in survival analysis or perhaps the third or fourth year undergrad stats course. This does appear to be a business-related task and maybe you should be thinking about soliciting a consultant's time? – DWin Aug 19 '16 at 23:30
• @DWin Thanks for the insight, it is indeed a business-related question, sadly what happened is that I'm currently in a internship where an actuary before taking a leave told me to do this myself, he also said that using the Cox model could do the work, but he didn't really look into the data himself it was just a guess. I find myself not able to complete the task I was given and I can't get any feedback, since he left until the end of my internship. I thought that maybe the situation was simple enough that someone here could explain to me how to solve it. I guess I just can't do this. – PaoloH Aug 22 '16 at 7:49

Well the coding answer IMHO is to use coxph on the Surv object, with the right hand side of the formula includes the factors and (maybe) other variables of interest. So something like ...

   coxph.fit <- coxph(Surv(time2013,event2013) ~ <contract type> + <other variables>, data=mydata)
summary(coxph.fit)


Where <...> indicates the actual factor/variable names.

For just two levels of a factor the resulting summary compares one level with the reference level. The p-value and confidence interval gives you maybe some idea of 'enough-ness' (i.e. statistical significance), but then you are getting into stats theory.

There's lots behind all this, so a bit of reading might be useful, to be aware of the assumptions and limitations.....

and maybe the guides here...

http://www.ats.ucla.edu/stat/r/examples/asa/asa_ch3_r.htm

Or follow a few threads here, to see the kinds of issues people have.

• I'll try to look into it again ,but I did try to use cox and last time what i got was that, even though there is a link between the number of sick leaves that are still going on after the x-th day and the number of days that happens after the x-th day, the results in the influence of variables did not seem to match. I'll post my results either case to let u know if it worked. – PaoloH Aug 18 '16 at 10:29
• I still have an issue with using Cox model so I edited my results in, and what makes me believe this might not be working for me. – PaoloH Aug 18 '16 at 14:18
• Some random observations in no particular order....1) the code formulation look okay to me, 2) you need to read up about proportional hazards. The plots don't show this - they cross over, 3) if you include CIs in the plots (conf.int), this would help support conclusions about differences, and should agree with results from summary(), 4) you also need to read up on Beta versus exp(Beta) if you haven't already done so, 5) you could plot with a log x-axis, 6) you might want to include censor ticks (mark.time) on the plot to check they aren't crazy. – Big Old Dave Aug 18 '16 at 15:47
• PaoloH - this is in response to @DWin's suggestion above (I cannot comment there due to low rep), but if you look at pages 14-15 of the first linked document above, it might get you started on parametric analysis using AFTs. If I was in your situation I would work through those examples in DIez, then apply your data when you feel more confident. It seems to me that you're also duty bound to report on the results obtained using CoxPH (as you've been asked/told to), but your conclusions should point out that the PH assumptions are not met. – Big Old Dave Aug 22 '16 at 18:44