5
$\begingroup$

I used the survival package in R to calculate the maximum likelihood estimates of Weibull regression co-efficients in R. I then tried to replicate the results by writing the log likelihood function and maximizing it using the maxLik package, using the "NR" method i.e Newton-Raphson. The point estimates from the two packages are similar, but the standard errors are extremely different (see below). I tried to investigate by finding the Hessian matrix for both approaches and found that they are also very similar, but not EXACTLY the same (see below). Am I doing something wrong in arriving at this conclusion that the Hessian is extremely sensitive? If I am right, then how does one get similar results computationally? If I want to maximize likelihood and find ML estimates of a custom function which is not available in any package, then how do I trust the standard errors that I get if the Hessian is so sensitive?

survival package point estimates: survival package point estimates

maxLik package point estimates:
maxLik package point estimates

Covariance matrix from both packages are very different, as shown below.

survival package covariance matrix: enter image description here

maxLik package covariance matrix: enter image description here

I found the Hessian for the survival package results by taking the negative of the information matrix, i.e the negative of the inverse of the covariance matrix). The Hessian matrices are almost the same, as shown below.

survival package Hessianenter image description here maxLik package Hessianenter image description here

EDIT FYI. Here is the output of vcov(fitWeibull). enter image description here

2nd EDIT adding information asked in the comment below

Basically, I am trying to reproduce the results from an example problem in a statistics text book Statistical Methods for Reliability Data, Second Edition (William Q. Meeker, Luis A. Escobar, Francis G. Pascual). I was able to do this using the survival package in R. I wanted to write the log-likelihood function and maximize it using the maxLik package. Again, I was able to reproduce the same point estimates. But the variance for the estimates that I get from the maxLik package is very different. The data consists of 26 observations (22 exact and 4 right censored) from a laboratory test, recording the number of cycles to failures for a given level of stress. The dataset is publicly available here. Since a curvilinear relation was observed in log(Thousands of Cycles) vs log(Stress), a quadratic term (square of log(Stress)) was added. Given below are text outputs that can be copied

survival package output

Call:
survival::survreg(formula = Surv(Cycles, delta) ~ logStress + 
    logStress2, data = superAlloy, dist = "weibull")
              Value Std. Error     z       p
(Intercept) 217.611     62.132  3.50 0.00046
logStress   -85.522     26.546 -3.22 0.00127
logStress2    8.483      2.831  3.00 0.00273
Log(scale)   -0.982      0.179 -5.48 4.2e-08

Scale= 0.375 

Weibull distribution
Loglik(model)= -93.4   Loglik(intercept only)= -118.4
    Chisq= 50.01 on 2 degrees of freedom, p= 1.4e-11 
Number of Newton-Raphson Iterations: 8 
n= 26

maxLik package output

m1 <- maxLik(loglik1, start = c(beta.0 = 0,
+                               beta.1 = 0,
+                               beta.2 = 0,
+                               sigma = 1),
+             method = "NR")
There were 42 warnings (use warnings() to see them)
> 
> summary(m1)
--------------------------------------------
Maximum Likelihood estimation
Newton-Raphson maximisation, 27 iterations
Return code 8: successive function values within relative tolerance limit (reltol)
Log-Likelihood: -93.38188 
4  free parameters
Estimates:
        Estimate Std. error t value  Pr(> t)    
beta.0 217.50945    4.96291   43.83  < 2e-16 ***
beta.1 -85.47893    2.12971  -40.14  < 2e-16 ***
beta.2   8.47809    0.23668   35.82  < 2e-16 ***
sigma    0.37474    0.06704    5.59 2.28e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
--------------------------------------------

I am confident that the log-likelihood function (loglik1) in the above code is correct, since I get the same point estimates. Notice the similarities in the point estimates and the huge difference in the standard errors above. I am trying to investigate this difference in the standard errors of the estimates. Here are the results of vcov for both the approaches:
survival package vcov

> vcov(fitWeibull)
              (Intercept)     logStress   logStress2  Log(scale)
(Intercept)  3860.3720980 -1649.1664686 175.82401657  0.29933426
logStress   -1649.1664686   704.7037272 -75.14987466 -0.12281312
logStress2    175.8240166   -75.1498747   8.01604061  0.01243889
Log(scale)      0.2993343    -0.1228131   0.01243889  0.03204363

maxLik package vcov

> vcov(m1)
              beta.0        beta.1        beta.2         sigma
beta.0  2.463048e+01 -10.368458628  1.0888625112  0.0007074572
beta.1 -1.036846e+01   4.535647200 -0.4950344210  0.0015390151
beta.2  1.088863e+00  -0.495034421  0.0560154009 -0.0004033332
sigma   7.074572e-04   0.001539015 -0.0004033332  0.0044947564

Since the survival package vcov consists of log(scale) instead of scale, I need to convert this matrix. I get the var(sigma) by multiplying the var(log(sigma)) by the square of sigma. I get the covariance(sigma, other point estimate) by multiplying the covariance(log(sigma), other point estimate) by sigma. This is by utilizing the delta method, which allows us to calculate the covariance of a function of random variables. Here is the result, which matches the results in the textbook:

> varCov.fitweibull
              (Intercept)     logStress    logStress2        sigma
(Intercept)  3860.3720980 -1.649166e+03 175.824016565  0.112172477
logStress   -1649.1664686  7.047037e+02 -75.149874664 -0.046022972
logStress2    175.8240166 -7.514987e+01   8.016040613  0.004661347
sigma           0.1121725 -4.602297e-02   0.004661347  0.004499885

I tried to investigate by comparing the information matrix (negative of the Hessian) from the 2 approaches, and here is what I get.

> solve(varCov.fitweibull) #information matrix for survival package "fitWeibull"
            (Intercept)  logStress logStress2     sigma
(Intercept)   156.66173   726.9912  3379.2575  29.61647
logStress     726.99120  3379.2575 15734.4024 140.35098
logStress2   3379.25748 15734.4024 73388.0562 666.03976
sigma          29.61647   140.3510   666.0398 229.46671
> solve(vcov(m1)) #information matrix for maxLik package "m1"
           beta.0     beta.1     beta.2     sigma
beta.0  156.68888   726.9989  3379.2276  29.64384
beta.1  726.99891  3379.3555 15734.1020 140.36061
beta.2 3379.22756 15734.1020 73384.7543 665.84960
sigma    29.64384   140.3606   665.8496 229.50530

As you can see, the information matrix is almost similar, but inverting them to get the covariance matrix gives totally different results. My question: "Am I doing something wrong in arriving at this conclusion that the Hessian is extremely sensitive? Or is my investigation wrong? If I am right, then how does one get similar results computationally? If I want to maximize likelihood and find ML estimates of a custom function which is not available in any package, then how do I trust the standard errors that I get if the Hessian is so sensitive?"

Improve this question
$\endgroup$
4
  • 2
    $\begingroup$ Try vcov(fitweibull) directly. Given the near density of the Hessian matrices, I suspect that the difference in the displays has to do with a $\sqrt N$ factor. $\endgroup$
    – EdM
    Commented Oct 26 at 17:29
  • $\begingroup$ I did use vcov(fitweibull). The default output of vcov(fitweibull) gives the covariance matrix with log(sigma). Since I wanted the covariance matrix with sigma, I converted the default matrix covariances to the one above using the Delta method. If you notice, the variance of sigma in both approaches actually match. It is all the other variances and covariances that are different. $\endgroup$
    – shishir rao
    Commented Oct 26 at 21:03
  • $\begingroup$ It would help if you could edit the question to show your results in text that can be copied, instead of images. You can do that by enclosing each block of results between two lines that start with three backticks. Also, I'd suggest that you say more about the data, your model, and how well it fits. You certainly seem to have some ill-conditioned matrices, and I wonder if that has something to do with the model itself. E.g., are your two predictors something like $x$ and $x^2$ instead of orthogonal polynomials? Details on how you "converted the default matrix covariances" would also help. $\endgroup$
    – EdM
    Commented Oct 27 at 17:52
  • $\begingroup$ Thank you @EdM. I have edited my question above and have added the information you requested. $\endgroup$
    – shishir rao
    Commented Oct 28 at 19:38

2 Answers 2

Reset to default
3
$\begingroup$

You are correct that the Hessian matrices and the covariance matrices derived from them for your models are highly "sensitive." What seems to work, at least in this example, is to improve the way the predictors are presented to the model.

The data are in the "SuperAlloy.csv" file at the linked data site. I was able to reproduce your original Weibull fit with survreg() after I multiplied the "Thousands of Cycles" values by 1000 and added extra columns containing log(Stress) and its square (not shown).

The variance-covariance matrix for that model has a large condition number, meaning that it is very close to a singular matrix.

kappa(vcov(fitWeibull),exact=TRUE)
# [1] 351742513

That can occur when there is severe multicollinearity among predictors. In this case, you can get around the problem by using an orthogonal polynomial in logStress instead of the separate terms:

fwOrthog <- survreg(Surv(Cycles,Fail)~poly(logStress,2),data=saData)
summary(fwOrthog)
# Call:
# survreg(formula = Surv(Cycles, Fail) ~ poly(logStress, 2), data = saData)
#                       Value Std. Error      z       p
# (Intercept)         10.8280     0.0821 131.91 < 2e-16
# poly(logStress, 2)1 -6.1099     0.4223 -14.47 < 2e-16
# poly(logStress, 2)2  1.3134     0.4384   3.00  0.0027
# Log(scale)          -0.9815     0.1790  -5.48 4.2e-08
# 
# Scale= 0.375 
# 
# Weibull distribution
# Loglik(model)= -245.4   Loglik(intercept only)= -270.4
#   Chisq= 50.01 on 2 degrees of freedom, p= 1.4e-11 
# Number of Newton-Raphson Iterations: 8 
# n= 26 

kappa(vcov(fwOrthog),exact=TRUE)
# [1] 37.55083

I'm not sure exactly how you handled the transformation of the covariance matrices to work with scale instead of log(scale), but I suspect that if you start with this much better conditioned matrix you will end up with much less "sensitive" results when you do that.

This page and its links discuss the advantages of orthogonal polynomials.

Improve this answer
$\endgroup$
1
0
$\begingroup$

Since all the point estimates were similar, I manually calculated the information matrix, taking partial derivatives and partial second derivatives of the log likelihood and estimating their values at the point estimates of the parameters. The covariance matrix that I got by inverting the manually calculated information matrix matched the results from the survival package and the textbook. Now, I had 3 different sources/methods where the covariance matrix matched (textbook, survival package and my manual calculation taking partial derivative) and 1 source/method where the covariance matrix did not match (maxLik package). I was a little uncomfortable putting this down solely to the highly sensitive Hessian matrix, so I took some time to read the maxLik package documentation. It strongly advises that users supply gradients (score functions/partial derivatives of the log-likelihood with respect to the parameters being estimated). Snip from maxLik package documentation

Once I added the gradient to the maxLik function, I was able to get the same standard errors as I get from all other methods. Here is the output of the maxLik with and without gradient:

maxLik without gradient

> loglik1 <- function(theta){
+   beta.0 <- theta[1]
+   beta.1 <- theta[2]
+   beta.2 <- theta[3]
+   sigma <- theta[4]
+  
+  logl <- sum(((log(failures.df$Cycles) - beta.0 - (beta.1*failures.df$logStress) - (beta.2*failures.df$logStress2))/sigma)) -   sum(exp((log(failures.df$Cycles) - beta.0 - (beta.1*failures.df$logStress) - (beta.2*failures.df$logStress2))/sigma)) - sum(log(sigma*failures.df$Cycles)) - sum(exp((log(censored.df$Cycles) - beta.0 - (beta.1*censored.df$logStress) - (beta.2*censored.df$logStress2))/sigma))
+     
+  return(logl)
+ }
> 
> m1 <- maxLik(loglik1, start = c(beta.0 = 0,
+                               beta.1 = 0,
+                               beta.2 = 0,
+                               sigma = 1),
+             method = "NR")
There were 42 warnings (use warnings() to see them)
> 
> summary(m1)
--------------------------------------------
Maximum Likelihood estimation
Newton-Raphson maximisation, 27 iterations
Return code 8: successive function values within relative tolerance limit (reltol)
Log-Likelihood: -93.38188 
4  free parameters
Estimates:
        Estimate Std. error t value  Pr(> t)    
beta.0 217.50945    4.96291   43.83  < 2e-16 ***
beta.1 -85.47893    2.12971  -40.14  < 2e-16 ***
beta.2   8.47809    0.23668   35.82  < 2e-16 ***
sigma    0.37474    0.06704    5.59 2.28e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
--------------------------------------------
> vcov(m1)
              beta.0        beta.1        beta.2         sigma
beta.0  2.463048e+01 -10.368458628  1.0888625112  0.0007074572
beta.1 -1.036846e+01   4.535647200 -0.4950344210  0.0015390151
beta.2  1.088863e+00  -0.495034421  0.0560154009 -0.0004033332
sigma   7.074572e-04   0.001539015 -0.0004033332  0.0044947564

maxLik with gradient

> loglik1 <- function(theta){
+   beta.0 <- theta[1]
+   beta.1 <- theta[2]
+   beta.2 <- theta[3]
+   sigma <- theta[4]
+  
+   z.F <- (log(failures.df$Cycles) - beta.0 - (beta.1*failures.df$logStress) - (beta.2*failures.df$logStress2))/sigma
+   
+   z.C <- (log(censored.df$Cycles) - beta.0 - (beta.1*censored.df$logStress) - (beta.2*censored.df$logStress2))/sigma
+   
+   #gradient for beta.0
+   grad.beta.0 <- sum((-1/sigma) - ((-1/sigma)*(exp(z.F)))) - sum(((-1/sigma)*(exp(z.C))))
+   
+   #gradient for beta.1
+   grad.beta.1 <- sum((-failures.df$logStress/sigma) - ((-failures.df$logStress/sigma)*(exp(z.F)))) - sum(((-censored.df$logStress/sigma)*(exp(z.C))))
+   
+   #gradient for beta.2
+   grad.beta.2 <- sum((-failures.df$logStress2/sigma) - ((-failures.df$logStress2/sigma)*(exp(z.F)))) - sum(((-censored.df$logStress2/sigma)*(exp(z.C))))
+   
+   #gradient for sigma
+   grad.beta.sigma <- (-1/sigma)*sum((z.F) - (z.F*exp(z.F)) + 1) + (1/sigma)*sum((z.C*(exp(z.C))))
+   
+   grad <- cbind(grad.beta.0, grad.beta.1, grad.beta.2, grad.beta.sigma)
+   
+   ll <- sum(z.F - exp(z.F) - log(sigma*failures.df$Cycles)) - sum(exp(z.C))
+   
+   attr(ll, 'gradient') <- grad
+   
+   ll
+     
+ }
> 
> m1 <- maxLik(loglik1, start = c(beta.0 = 0,
+                               beta.1 = 0,
+                               beta.2 = 0,
+                               sigma = 1),
+             method = "NR")
There were 34 warnings (use warnings() to see them)
> 
> summary(m1)
--------------------------------------------
Maximum Likelihood estimation
Newton-Raphson maximisation, 13 iterations
Return code 1: gradient close to zero (gradtol)
Log-Likelihood: -93.38188 
4  free parameters
Estimates:
        Estimate Std. error t value  Pr(> t)    
beta.0 217.61114   62.13232   3.502 0.000461 ***
beta.1 -85.52238   26.54644  -3.222 0.001275 ** 
beta.2   8.48273    2.83128   2.996 0.002735 ** 
sigma    0.37474    0.06708   5.586 2.32e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
--------------------------------------------
> vcov(m1)
              beta.0        beta.1        beta.2        sigma
beta.0  3860.4254598 -1.649189e+03 175.826442951  0.112177097
beta.1 -1649.1892470  7.047135e+02 -75.150910408 -0.046024946
beta.2   175.8264430 -7.515091e+01   8.016150941  0.004661557
sigma      0.1121771 -4.602495e-02   0.004661557  0.004499886

The covariance matrices now match!

The fact that the point estimates from the maxLik package without providing the gradient matched other sources combined with the fact that the Hessian was also very similar to other sources led to my line of investigation focusing on the Hessian. It is also true that the Hessian matrix can indeed be very sensitive, as EdM has helpfully pointed out in the other answer.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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