2011:

For the Empirical Likelihood ratio (Wilks) theorem for hazard integrals with
censored data: show the two modifications:

(a) Verify that the null hypothesis  \int g(t) d H(t) = \theta_0, can be replaced by
    \int g(t) d H(t) = \theta_0 + o(1/sqrt n) and the chi square
   limit for the -2 likelihood ratio still hold under H_0.

(b) If we replace the null hypothesis by \int g(t) dH(t) = \theta_0 + ( C / \sqrt n)
    then the limit of the -2 log likelihood ratio has limit of a non central 
    chi square distribution. Identify this non-centrality parameter.

You may want to consult the note 2.1. TO make things simpler, you may assume just
one constraint. (So the resulting chi square has df=1)

==============================================================================

Is there a simple proof for the NPMLE of Kaplan-Meier? like Owen 2001's
proof of the sampling distribution being NPMLE?


Added 2010 Spring:

Exercise: Without extra constraint, the maximization of the nonparametric likelihood
function with right censored data, is solved by Kaplan-Meier (1958). Try to solve
it yourself, without consult the paper.  
There is an explicit solution given by Kaplan-Meier,
there is also a recursive solution. 



For 2007 Spring

Jan.

(0) Use a data set with more than 100 observations.
Use SAS proc lifetest (version 9)
obtain the confidence band for the Kaplan-Meier
curve and point-wise confidence interval. Compare.

Try to move the upper and lower time boundary for the 
confidence band closer to each other, observe the changes.

Some keywords for SAS: confband = ep hw all
                       bandmax
                       bandmin
                       conftype 

For people using R, you can do the same thing with
a package called km.ci  Load the package and read the help
file. 

(1) Suppose BM(t) is a Brownian Motion process.
Define the process
BB(t) = BM(t) - t BM(1),  for 0 <= t <= 1

Compute the variance of BB(t), and compute
the covariance between BB(t) and BB(s) for
0 <= t < s <= 1.

Verify that BB(t) do not have independent increment.

=============================================================
(2). Suppose X_1, X_2, ... X_n are iid r.v.s with a cont. CDF F(t).

Let  U_n(t) = \sqrt n [ \hat F_n(t) - F(t) ]

Compute the covariance of U(t) and U(s), i.e. E U(t) U(s) = ?
(where \hat F_n(t) is the empirical distribution based on X_i's).

(3). If N(t) is a Poisson process, show that N(3t) is also a Poisson
process.



(4) Suppose N_i(t), i=1,2, ... , n are iid Poisson ( lambda )
processes.   0 < t < 8 . (or some other fixed finite number).

Show that: \sum N_i(t) is also a Poisson process with parameter (n lambda)

So \sum N_i(t) and N_1( n t) are both Poisson (n lambda) process.

=====================================================================

Feb. 

(5) Let M(t) = N (t) - t  where N(t) is a standard Poisson process
and also let us define the filtration
F_t = history sigma algebra for info in [0, t]

We verified in class that M(t) is an F_t martingale.

Please verify  [M(t)]^2 - t   is also an F_t martingale.
(hint, the variance of Poisson r.v. is also lambda)


March.   

Use the technique we learned in class to derive
the asymptotic variance of 

\sqrt n  [ \int_0^\infty  g(t) d \hat F(t) - \int_0^\infty g(t) dF(t)] 

assume the said variance is finite. Here g(t) is a continuous differentiable
function.

Hint: first use integration by parts to turn the above into

\sqrt n [ \int_0^\infty 1-\hat F(t) dg(t) - \int_0^\infty 1-F(t) dg(t)]

here \hat F(t) is the Kaplan-Meier estimator.



A related but different quantity: find the asymptotic
variance of the follow

\sqrt n [ \int_0^\infty g(t) (1-\hat F(t-))d \hat H(t) - 
          \int_0^\infty g(t) (1-\hat F(t-)) dH(t) ]

here H(t) is the cumulative hazard function, and \hat H(t) is the
Nelson-Aalen estimator.

April

Make use of the asymptotic variance formula we derived in March, of the quantity
\sqrt n [ \int_0^\infty 1-\hat F(t) dg(t) - \int_0^\infty 1-F(t) dg(t)],

derive an estimator of the asymptotic variance.
 
Now construct a test based on the mean of g(): standardized to 
have variance 1

simulate the distribution (histogram) of the normalized
test statistics and confirm it is approx. N(0,1)
for samples of size 60 and (around) 15% censoring. use the following
g() functions:  g(t) = exp( -0.1 t^2 )
                g(t) = t for   0<= t <= 3
                     = 3 for       t > 3

You may want to take a look of my preliminary R code  
www.ms.uky.edu/~mai/sta709/varG.R
for estimating the variance.


Remark: if you do not get exactly 15% censoring, do not worry.
anything between 10% and 25% is OK.

I realize there are more problems with the coding of this ..... please take a
look of the code in the above mentioned place. I updated the code.



Last homework: A two sample test, based on the Kaplan-Meier mean: 
In the previous problem we worked on the one sample. 
If there are two samples, then a similar "t-test" can also be worked out.
Testing whether the two (censored) samples have the same mean.
Variance would be the sum of the two (by independence).
i.e.

numerator:     \int g(t) d \hat F_1(t) - \int g(t) d\hat F_2(t)

denominator:   \sqrt ( Var est. 1 + Var est. 2 ).

Compare this with the (2-sample) log-rank test.




Interval censor: some code 
(Reference: The Statistical Analysis of Interval-Censored 
Failure Time Data, published by Springer, 2006 author: J. Sun)
page. 77 


simuIC <- function(n) {
x <- 1:n/n                     ## x is the observation times, 
xorder <- order(x)             ## ordered 
x <- x[xorder]                 ##
z <- rbinom(n, size=1, p=0.5)  ## z is the indicator of the sample1 or 2
z <- z[xorder]                 ## z also ordered according to x
y1 <- rexp(n)                  ##
y2 <- rexp(n)                  ## y1 and y2 are two samples of lifetimes
y <- y1*z + y2*(1-z)           ## put the two samples together
yy <- as.numeric( y < x )      ## interval censoring, current status
temp <- isoreg(x=x, y=yy)      ## get the estimator hat F, as yf in temp
weights <- runif(n)            ## create some random weights
SDF <- sum(z*(yy - temp$yf)*weights)      ## this is the numerator of the test
varSDF <- var(z*weights)*sum((yy - temp$yf)^2)  ## var est of the numerator
ZZ <- SDF/sqrt(varSDF)         ##  test statistic
list( x=x, y=yy, f= temp$yf, ZZ=ZZ, diff=SDF, VD=varSDF)   ## output list
}






Takehome final problem:

A quantity important in some analysis is called AUC 
(Area Under Curve, for ROC curves). 

The definition of the AUC goes like this: 
There are two independent samples of survival times, with
sample size n and m; and  
with population survival distributions F_1 and F_2 respectively. 

AUC = \int F_1(t) d F_2(t) . 


Based on the right censored data/sample,
we can get two Kaplan-Meier estimators for the two samples.
An estimator of AUC is obviously to replace the unknown F by the hat, 
(For censored data, \hat F should be the Kaplan-Meier estimator)

\hat AUC = \int \hat F_1(t) d \hat F_2(t).

What is the asymptotic variance of (as in a formula)

\sqrt(n+m) [ \hat AUC - AUC ] ?

assume both n and m increase to infinite, and the 
ratio n/m remains not close to 0 and not close to 1. 

Assume the variance is finite.


Hint:  (1) show that \hat AUC - AUC = 
           \int (\hat F_1 - F_1) dF_2 + \int F_1 d(\hat F_2 - F_2) + \epsilon
           (argue, at least heuristically \epsilon is the high order 
            small term, involve at least two differences, or 
            difference square...)

       (2) Assume \sqrt(n+m) \epsilon is negligible as sample sizes grows,
           we focus on the two integrals. Variance of first integral can 
           be handled similar to homework.
           Integration by parts on the second integral will bring it to 
           the familiar form. (assume you always define the last obs. as 
           uncensored)


Anyone interested on some further reading, please see the comments
at the bottom of this file.



















Old homework
==============================================

(1). Suppose N(t) is a Poisson process.

Verify the following is a (continuous time) martingale.

M(t) = N(t) - lambda t 

for t in [0, T]

Be sure to define your filtration (increasing family
of sigma algebras) F_t. (The history)

====================================================

(2). Suppose N_i(t), i=1,2, ... , n are iid Poisson ( lambda )
processes.   0 < t < 8 . (or some other fixed finite number).

Show that

     M_n(t) =  \sum [ N_i(t) - lambda t ] / sqrt (n lambda )

as n \to \infty has the following property

    (a) for fixed t the distribution of M_n(t)
        converges to a Normal distribution with 0 mean
        Identify the variance.

    (b) the sample path of M_n(t) has jumps of size converge to 0. 
           (1/sqrt n \lambda)

   Note: \sum N_i(t) is also a Poisson process with parameter (n lambda)

=======================================================
(3). If N(t) is a Poisson process, show that N(3t) is also a Poisson 
process.
===============================================================
(4). Suppose X_1, X_2, ... X_n are iid r.v.s with a cont. CDF F(t).

Let  U_n(t) = \sqrt n [ \hat F_n(t) - F(t) ]

Compute the covariance of U(t) and U(s), i.e. E U(t) U(s) = ?
(where \hat F_n(t) is the empirical distribution based on X_i's).

================================================================
(5).  Suppose M(t) is the Poisson martingale we defined in the first problem.
Show that 

      M^2(t) - \lambda t  

is again a (cont. time) F_t martingale. 
(same filtration we used there). 

(5.5) Suppose M(t) is a martingale. Show that
       E[ M^2 (t + dt) - M^2 (t) |F_t] = E[ ( M(t+dt) - M(t) )^2  |F_t] 
====================================================================
(6) Let N(t)= I[ X <= t ], where X is a positive, cont. rv.
    Define the filtration F_t = \sigma { N(s); 0 <= s <= t }.
    (a) show that  h(t)I[X >= t] is predictable.
        where h(t) is the hazard function of X.

    (b) verify the conditional expectation of the 
            increment of N(t) over (t-dt, t] , given F_{t-dt} 
            (i.e. either given I[X >= t]=0 or given I[X >= t]=1)
            is h(t)I[X >= t] dt 

    Note this shows that N(t) - \int_0^t h(s)I[X >= s] ds
        is an F_t martingale.
==================================================================
(7) We proved the martingale representation of the Kaplan-Meier estimator
for continuous CDF F(). Now Suppose that the CDF F() hava a jump at t, 
verify the
martingale representation of the Kaplan-Meier estimator is still valid at t.

For simplicity, assume all the observed death times do not equal to t.
Also, when the CDF F() is not continuous, the compensator
of the counting preocess needs to be written as

\int_0^t I_{[X >= s]} d H(s) ;  instead of

\int_0^t I_{[X >= s]} h(s) ds .
=================================================================
(8)
Consider a stochastic jump process N(t), t>=0 that fit the 
following description:
(a) N(0) = 0 
(b) it can have at most one jump.
(c) N(t) have probability zero to have jump for 0 <= t <= a
(d) after the point a, the waiting time to the jump is 
Weibul( 1, 0.2) distribution;
(e) the size of the jump is inversely proportional to the 
waiting time after a.
(f) if after time c, N(t) still have no jump, then it will never jump.

Assume 0<a<c are known constants. Write down the 
compensator/cumulative intensity of N(t)
so that N(t)- A(t) is a martingale wrt 
$ {\cal F}_t= \sigma \{ N(s), 0 <= s <= t \}$

===========================================================================
(9) Using the martingale representation for the Kaplan-Meier estimator
\[
      \frac{\hat F_n (t) - F(t)}{1-F(t)} = 
     \int_0^t \frac{1-\hat F_n(s-)}{1-F(s)} \frac{1}{Y(s)} d \sum M_i(s)
\]
  
to find out the covariance
\[
E n[\hat F_n(t) - F(t)][ \hat F_n(s)-F(s)]
\]
simplify the expression as $n \to \infty$.
============================================================================
(10) Suppose $\hat F_n(t)$ is the empirical CDF based on n iid observations.
Explain why $\hat F_n(t) - F(t)$ cannot have a martingale representation:
$ \int_0^t H(s) dM(s)$.
Hint: consider variance and predictable variation process.
=============================================================================
(11) 

Consider a 
pure jump stochastic process N(t), for t >= 0  (sample path 
are piecewise constant)
that fit the following description: 
(a) N(0) == 0
(b) it can have at most two jumps.
(c) N(t) have probability zero to jump for 0 <= t <= a
(d) after the point a, the waiting time to the first jump is
Weibul(2, 0.2) distribution; The waiting time to the second jump 
(after the first jump occured) is exp( \lambda= length of first waiting time). 
(e) after time c, N(t) will never jump.

Assume 0<a<c are known constants. Write down the
compensator/cumulative intensity of N(t)
so that N(t)- A(t) is a martingale wrt
$ {\cal F}_t= \sigma \{ N(s), 0 <= s <= t \}$

(f) Everything else the same as before but the 
stochastic jump process have size of the 
first jump inversely proportional to the waiting time of first jump,
the size of the second jump equal to a random variable Y. Assume Y 
independent of everything else.
Use $ {\cal F}_t= \sigma \{ Y, N(s), 0 <= s <= t \}$ here.
============================================================
(12) 
This question refers to the data set at 

http://www.mayo.edu/hsr/people/therneau/book/data/gastric.html


1). There are two treatments. For each treatment, 
at t= 365, use empirical likelihood
method [ function emplikH1.test() in R package emplik ], 
to obtain two confidence intervals
of H_i(365) where H_i() is the cumulative hazard function of trt i

Does the confidence intervals overlap?

2). Now use the function g(t) = exp(- 0.003*t) and get two confidence
intervals of 

     \int g(t) dH_i(t) 

by empirical likelihood method.

Does the confidence intervals overlap?

3). Use Wald method to obtain confidence intervals of H_i(365).
and compare to those obtained in 1). 

All the confidence level = 95%.
======================================================
Answer to the last question:

The 95% confidence interval for H1(365) should be

    [ 0.2011, 0.6068 ]  by emplik method. 

The 95% confidence interval for H2(365) should be 

    [ 0.5194, 1.1657 ]  by emplik method.

The H2(365) is higher but not significant.

The 95% confidence interval for   int_0^\infty  exp(-0.003*t) dH1(t)
should be 
             [ 0.2715, 0.5493 ]  by emplik method.

The 95% confidence interval for   int_0^\infty  exp(-0.003*t) dH2(t)
should be 
             [ 0.4505, 0.8743 ]  by emplik method.

Again int_0^\infty  exp(-0.003*t) dH2(t)  is higher but not significant.


Use Wald method, we get 95% confidence interval for H1(365) as

             [ 0.1713, 0.5640 ]

and the 95% confidence interval for H2(365) as

             [ 0.4705, 1.1239 ] . 
===========================================================
If M(t) is a counting process martingale wrt F_t.
Prove or give a counter example that

(a)  g(t)M(t) is also a martingale. where g(t) is an increasing function.
(b)  3M(t) is also a martingale.
(c)  3N(t) is a Poisson process if N(t) is. 
=============================================
For Poisson process, the inter-arrival times are iid exponential r.v.s.
Suppose a counting process whose inter-arrival times are iid
U(0,1) r.v.s

Find the intensity function of this counting process.
===============================================
How do you Simulate
a counting process with intensity function sqrt t ?
===============================================
How do you simulate a pure jump stochastic process
(one that is constant between jumps): the time between jumps 
are iid exponential (1), but the size of the jumps
are: 1 for the first jump, 1/2 for the second jump,
1/3 for the third jump, ... 1/n for the nth jump.

Find also the compensator for this process so that
the difference is a martingale.  
====================================================
what if in the above the jump sizes are exp(-t) where
t is the time of the jump.
========================================================================
what if the interarrival time is exp(k) for the time
between k and k+1 th jump?
Use simulation to find/compare the performance
of the two/three types of censored mean EL confidence interval.

===============================================
Use simulation to compare the two/three types of logrank test.
============================================================
Suppose Y_i are independent exp( \lambda_i) r.v.s, i=1,2,...n

Find the compensator for the counting process

N(t) = \sum I [ sqrt(Y_i) \leq t ] 

so that N(t) - A(t) becomes a martingale.
===============================================
Suppose a pure jump stochastic process N(t) has interarrival time
distributed as exp(k) for the time between k and k+1 th jump
(and independent).  The size of the jump is 1/k. 
[The first jump has size 1 and waiting time exp(1).]

Find the compensator A(t) of this stochastic process so that
N(t) - A(t) is a martingale.
============================================================
The power of the logrank type test:

(1) Two sample testing case. Sample one from an
    exp(lambda) distribution. Sample two from a distribution
    with hazard(t) = B exp( - a t )  (cross hazard alternative.)

    Show by simulation example that for certain choice of B and a, 
    the log rank test have bad power for this test.
    Use sample size 200 and 200. 
    (survdiff, phreg, and emplikHs.test2 ) 

(2) the logrank test statistic can be thought of as an
    integration over all the time range.
    How can you derive a test that only integrate
    up to a fixed time T? (how to do it in R or SAS?)
    Integraltion over an interval [U, T] ?

     Does the power improve if we only compute the logrank test
     over certain interval for data as in (1)?

=================================================================

Some commemts related to the takehome final problem:
Background reading:  http://www.bepress.com/uwbiostat/paper268
(just the first half of the paper is enough).
 

How to maximize the 2-sample Empirical Likelihood
under the constraint of AUC = 0.8 (say)?
Assume no censoring for the moment.
So the observations are x1, x2, ... xn from sample 1, and 
y1, y2, ... ym from sample 2. with population CDF F() and G().

The two sample Empirical Likelihood should be defined as

\prod_{i=1}^n p_i \prod_{j=1}^m w_j

The max of the EL without additional constraint is obviously
obtained when p_i=1/n and w_j = 1/m.


When maximize the probabilities under the constraint of AUC = 0.8, 
we write the probability as

     p_i = \frac{1}{ n + \lambda ( \hat G(x_i) - 0.8 ) }

and
     w_j = \frac{1}{ m - \lambda ( \hat F(y_j) - 0.2 ) }

Now try to iteratively solve the \lambda by (1) fix \hat F, you can
get w_j which is \hat G, and 
then (2) fix \hat G, you can get p_i, which is \hat F. etc. 
until convergence.  (May be you need to half the steps in search of lambda)
and got a solution \lambda. I think this lambda and the corresponding
p_i, w_j should achieve the max of the EL under constraint.
(needs verification)

Here we have one unknown lambda but two equations....the solution
is not straight forward.

Also, need to verify -2 log lik ratio has chi sq distribution under Ho

Assume both F and G are cont. so no ties in the data.


The next step is to allow both sample x and y to be right censored.
In this case, we do not have explicit lambda solutions for the probability
but we still have p_i (no constraint) and  p_i (with constraint).
The p_i (no constraint) is just the jumps of the Kaplan-Meier estimator
p_i (with constraint) can be obtained using the emplik package function
el.cen.EM.   In this case the constraint is  \int G(x) dF(x) =0.8
or E G(X) = 0.8.  Which G to use is a problem and has to be chosen as
the Kaplan-Meier based on the y_j and later updated in iteration. 

[With the new version of emplik 0.9-3-2, the function el.cen.EM will
also output a lambda value, this may be used to solve maximization problem]
Since lambda can be viewed as a "degree of deviation from NPMLE"
and in two samples, we want to have the two samples to have same degree
of deviation, one move to the right, one move to the left, so they 
should have the same lambda.


I think we shall not take p_i (with constraint) as is. Because
I think it try to correct the probability by adjust p_i only 
and thus is over-correct (without changing w_j).
we probably should take p^* = p + 0.5( p(with constraint) - p)
i.e. only took half the step. In any case, when convergence is achieved
we shall have p^* and w^* such that : use p^* or F^* we shall get w^*
from the el.cen.EM. and using w^* or G^* we should get p^*.
That means convergence. (this convergence may not be unique?)
[So use lambda value to test convergence]


Another possible generalization of Owen, (2001)
2 sample ELT, is to view the function
I[X>Y] as a kernal of a U-statistics. 
And make the ELT theorem hold for all U-statistics, 
(well, all those with order 2 and finite variance) 
for example H(x,y) = xy.  Or is this trivially true
given Owen's outline proof?

What about 3 samples?

For censored data, we can try I[X>Y] too. What about U-stat?