2
$\begingroup$

I wanted to check how fast an average colleague of mine is able to complete a puzzle, so I ran some experiments. The problem is they weren't done in exactly 50 minutes, I always got bored and moved on to the next person. Hence my data is right-censored at 50. Let's say it looks like this:

Aaron: 34
Betty: X
Cecil: X
Dave: 18
Elias: 47
Fredric: 39
Greta: 37
Howard: X
Ingrid: X
Joy: 29
Kristine: 34
Lou: 32
Margret: X
Nigel: 40

An X means that they did not complete it in 50 minutes.

How do I use this data to find how fast on average person in my company is in completing the puzzle? To be more clear: How do I find the mean or the median of the above data with a bootstrap confidence interval? I don't want to just throw out the data where the person did not complete the task within 50 minutes, but I don't know how to incorporate those data points in my analysis.

$\endgroup$
3
  • $\begingroup$ The mean can abviously not be computed without further assumptions and the median should be no problem, as long as less then half of the population has no X. Just replaxe X with a very large number (if you are doing this in R I suggest Inf) and compute the median. Depending on the number of X you might also consider a trimmed mean. $\endgroup$
    – Bernhard
    Mar 13, 2018 at 10:14
  • $\begingroup$ @Bernhard. Thanks. It may happen that more than half of the population has 'X'. Does that mean that I cannot draw any conclusions? $\endgroup$
    – Sid
    Mar 13, 2018 at 11:32
  • $\begingroup$ Intuitively, without math, what do you expect the mean of $(3, 5, 9, X, X, X, X, X, X)$ to be? More or less then $20$? $\endgroup$
    – Bernhard
    Mar 13, 2018 at 15:09

2 Answers 2

1
$\begingroup$

You can use interval regression on a constant, since you have a mixture of point data, which can be expressed as $[t_1,t_2]$ where $t_1=t_2 < 50$, and right-censored data, where completion time is $[50,+\infty)$. In this case, the expected duration is 43 minutes, instead of the 40 you get from ignoring the right-censored folks. This model relies on the error being normally distributed with mean zero and homoskedastic variance. This is the sort of strong parametric assumption that Glen_b is rightly concerned about.

Here's the example in Stata:

. clear

. input str8 name str2 time

          name       time
  1. "Aaron" 34
  2. "Betty" X
  3. "Cecil" X
  4. "Dave" 18
  5. "Elias" 47
  6. "Fredric" 39
  7. "Greta" 37
  8. "Howard" X
  9. "Ingrid" X
 10. "Joy" 29
 11. "Kristine" 34
 12. "Lou" 32
 13. "Margret" X
 14. "Nigel" 40
 15. end

. 
. gen t1 = cond(missing(real(time)),50,real(time))

. gen t2 = cond(missing(real(time)),.,real(time))
(5 missing values generated)

. intreg t1 t2, nolog

Interval regression                             Number of obs     =         14
                                                   Uncensored     =          9
                                                   Left-censored  =          0
                                                   Right-censored =          5
                                                   Interval-cens. =          0

                                                LR chi2(0)        =       0.00
Log likelihood =  -40.91248                     Prob > chi2       =          .

------------------------------------------------------------------------------
             |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       _cons |   43.20936   4.073829    10.61   0.000      35.2248    51.19391
-------------+----------------------------------------------------------------
    /lnsigma |   2.634876   .2589092    10.18   0.000     2.127423    3.142329
-------------+----------------------------------------------------------------
       sigma |   13.94158   3.609604                      8.393212    23.15773
------------------------------------------------------------------------------

If you still want to boostrap, that can also be done:

. bs, reps(100): intreg t1 t2, nolog
(running intreg on estimation sample)

Bootstrap replications (100)
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5 
..................................................    50
..................................................   100

Interval regression                             Number of obs     =         14
                                                   Uncensored     =          9
                                                   Left-censored  =          0
                                                   Right-censored =          5
                                                   Interval-cens. =          0
                                                Replications      =        100

                                                Wald chi2(0)      =          .
Log likelihood =  -40.91248                     Prob > chi2       =          .

------------------------------------------------------------------------------
             |   Observed   Bootstrap                         Normal-based
             |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       _cons |   43.20936    4.94787     8.73   0.000     33.51171      52.907
-------------+----------------------------------------------------------------
    /lnsigma |   2.634876   .2534669    10.40   0.000      2.13809    3.131662
-------------+----------------------------------------------------------------
       sigma |   13.94158    3.53373                      8.483219    22.91203
------------------------------------------------------------------------------
$\endgroup$
2
  • $\begingroup$ Cool, thank you! I was going to experiment with this myself, but I saw that Strata is proprietary software. Do you have any pointers on how I can do this in Python or R? $\endgroup$
    – Sid
    Mar 14, 2018 at 12:15
  • $\begingroup$ Try the censReg package in R. This simple case is also called a Tobit model. $\endgroup$
    – dimitriy
    Mar 14, 2018 at 16:30
3
$\begingroup$

The bootstrap (alone) will not help you identify the mean in the presence of censoring, and will not be much use with the median given a dataset like the one you supplied; you may be able to get some kind of lower bound on the median in this case with a bootstrap, but you can do that nonparametrically without the need for bootstrapping.

To estimate the mean would require some kind of parametric assumption*.

If you observe fewer than half the times to completion then the median would have a similar problem.

This is a standard kind of problem in survival analysis; once a suitable parametric distributional model is chosen, estimation is relatively easy with canned routines. I've done this in R; it provides a number of distributions by default but it's possible to add further distributions.

First we create a dataset with time under observation and a column to denote the censoring status (first few rows shown below). Here I fit an intercept-only Weibull model to the times

  subject time completed
1   Aaron   34         1
2   Betty   50         0
3   Cecil   50         0
4    Dave   18         1
5   Elias   47         1
6 Fredric   39         1

(weib.fit <- survreg(Surv.puz~1,puzzle,dist="weibull"))
Call:
survreg(formula = Surv.puz ~ 1, data = puzzle, dist = "weibull")

Coefficients:
(Intercept) 
   3.879441 

Scale= 0.3010924 

Loglik(model)= -40.8   Loglik(intercept only)= -40.8
n= 14 

Here's the Kaplan-Meier survival curve (nonparametric), along with 95% CI, and the fit from the above model:

Plot of Kaplan Meier survival curve with 95% CI and fitted Weibull, showing fitted median and mean

The sample median is 39.5, the fitted median is 43.34, the fitted mean is 43.43. You can see from the nonparametric CI that we can get a lower bound on the median without a parametric assumption (where the horizontal green dots hit the vertical black dashes) but there's not a large enough proportion of complete times to get an upper bound on the median. The parametric fit can provide CI-type bounds on both the mean and median (though I haven't given any here).

You could try a variety of parametric assumptions if you don't want to be too reliant on any one assumption, but you're still reliant on whatever set of distributions you consider -- there's no way to avoid some assumption about the upper tail beyond what you observed.

* well it might be possible to proceed with some assumption that sufficiently regularizes the problem without being fully parametric (obtaining some set of bounds, presumably), though I am not sure right now what form that assumption might take.

Here's the code that I used. First read in and set up data:

puzzle <- read.table(stdin(),colClasses=c("character","numeric"),
             na.strings="X",sep=":",col.names=c("subject","time")) 
Aaron: 34
Betty: X
Cecil: X
Dave: 18
Elias: 47
Fredric: 39
Greta: 37
Howard: X
Ingrid: X
Joy: 29
Kristine: 34
Lou: 32
Margret: X
Nigel: 40

 puzzle$completed <- as.numeric(!is.na(puzzle$time))
 puzzle$time <- with(puzzle,ifelse(is.na(time),50,time))

Then load the survival library and create the required survival object, fit the Weibull and plot the curves etc

 require(survival)
 (Surv.puz <- with(puzzle,Surv(time,event=completed)))
 (weib.fit <- survreg(Surv.puz~1,puzzle,dist="weibull"))

 plot(survfit(Surv.puz~1),xlim=c(0,100))
 lambda <- exp(weib.fit$coefficients)
 k <- 1/weib.fit$scale
 f <- function(x) pweibull(x,k,lambda,lower.tail=FALSE)
 curve(f,0,100,col="blue",add=TRUE)
 weib.mean <- lambda*gamma(1+1/k)
 weib.median <- lambda*log(2)^(1/k)
 puz.median <- median(puzzle$time)

 abline(h=c(0,1),col="darkred",lty=3)
 segments(0,.5,max(weib.median,puz.median),.5,col="green4",lty=2)
 segments(weib.median,0,weib.median,.5,col="green4",lty=2)
 segments(puz.median,0,puz.median,.5,col="green4",lty=2)
 arrows(weib.mean,-.03,weib.mean,0,length=0.05,col="green4")
 c(data.median=puz.median,weib.median=weib.median,weib.mean=weib.mean)
$\endgroup$
5
  • $\begingroup$ +1 I... I just can't help it! The white background of graphs clashes really hard with the faint graph-paper background of the site... I am repeatedly compelled to edit graphs for transparency... there's probably a 12-step program for this kind of thing. ;) $\endgroup$
    – Alexis
    Mar 14, 2018 at 0:39
  • $\begingroup$ @Alexis Well, I'll leave it as it is for your sake, but personally I prefer the white background to having the page background showing through -- it's not enough to bother me, but I like the "clean paper" look of the white background on plots. On my monitor the difference is pretty subtle though. $\endgroup$
    – Glen_b
    Mar 14, 2018 at 5:15
  • $\begingroup$ How are you able to get the diagram and the data (sample median and fitted median) using the above R code? $\endgroup$
    – Sid
    Mar 22, 2018 at 3:03
  • $\begingroup$ I didn't show every bit of R code I used, but this can be worked out from the parameters estimates above using information in the help of the relevant functions. In terms of the parameterization of the Weibull given at Wikipedia, $\hat\lambda=\exp(\tt{Intercept})$ and $\hat{k}=1/\tt{Scale}$, from which the maximum likelihood estimates for the median and mean may be computed by substitution. You can draw a Weibull survivior function via a simple call to pweibull. The Kaplan-Meier 'curve' is drawn via plot.survfit. $\endgroup$
    – Glen_b
    Mar 22, 2018 at 3:51
  • $\begingroup$ ... I only kept the code a few days; I'd have to regenerate it from scratch. Actually, I see a small error in my copy of your data, which I will also fix. (now done. Apologies to Alexis; I'll try to fix up the transparency later.) $\endgroup$
    – Glen_b
    Mar 22, 2018 at 4:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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