[R] Tobit Estimation with Panel Data

[Terry Therneau]

  Adding true random effects to survreg is certainly on my list of useful 
additions, but one with no start date in sight.  That said, one can get an 
alternate solution with
     survreg(Surv(time, status) ~ x1 + x2 + frailty.gaussian(id, method='aic'))

  Justification: one can view a random effects model as a penalized model, that 
is, as
  	a. the addition of "factor(id)" - a coefficient b_i for each subject
  	b. a shrinkage penalty,  -.5*k* sum(b_i^2), is added to the log-lik, and
  we minimize the sum
  	c. the value of k is chosen to maximize an integrated likelihood, one    
  with the b's integrated out.  1/k is the variance of the random effect
  	
The above code uses the AIC to choose k.  You could also use a user-specified 
degrees of freedom.

 WARNING: Using "method='reml'" in the above won't work correctly.  The fact 
that no warning is given in this case is serious flaw in the survival package. 
(In my defense, the 'penalty function' code for coxph and survreg was designed 
to allow general user-written penalties; a side effect is that the penalty 
functions can't easily tell which routine is calling them.  Most, e.g., 
pspline() and ridge(), work for both coxph and survreg.  But frailty with either 
an 'ml' or 'reml' argument computes the appropriate Cox model integral, not the 
survreg one, and so gives nonsense answers when used with survreg.)

	Terry Therneau