Struggling to understand measures of influence in Cox PH models I am currently reading a text book that recommends assessing how much influence any individual subject may have on a fitted cox regression model. In particular it recommends using displacement of log-likelihood after removing subject i as a measure of influence. That is to say
$$
LD_i = 2 ( L(\hat \beta) - L( \hat \beta_{(i)} ) 
$$
where
$L(.)$ is the log-likelihood function
$\hat \beta$ are the estimates of our cox-model coefficients
$\hat \beta_{(i)}$ are the estimates of our cox-model coefficients after removing subject i 
The text then goes on to say that an approximation of this value was shown by Pettitt and Bin Daud (1989) to be 
$$
LD_i = r_{U_i} var(\hat \beta) r_{U_i}^T
$$
Where
$r_{U_i}$ is the p x 1 vector of score residuals for subject i. 
My understanding then begins to fall apart when I compare this approximation back to the true value of the likelihood displacement (see graph and reproducible code below). In particular what seems odd to me is that there is a lack of a linear relationship / positive correlation between the 2 values and that they appear to be on completely different scales.  Am I misunderstanding what these values are supposed to be ? 

R Code:
library(survival)
library(ggplot2)

### Poor mans imputation to stop errors associated with coxph just dropping rows
lung[ is.na(lung$inst),]$inst <-  16

MOD <- coxph( 
    Surv(time, status -1) ~  age + inst, 
    data = lung 
)

SCORE <- residuals(MOD , "score")

LD_APPROX <- sapply ( 
    1:nrow(lung), 
    function(i) t(SCORE[i,]) %*% MOD$var %*% SCORE[i,]
)

get_ld <- function(i , MOD){
    tlung <- lung[-i,]
    MODW <- coxph( 
        Surv(time, status -1) ~  age + inst , 
        data = tlung 
    )
    abs(2 * as.numeric(logLik(MOD) - logLik(MODW)))
}

dat <- data.frame(
    LD = sapply( 1:nrow(lung) , get_ld , MOD = MOD),
    LD_APPROX = LD_APPROX 
) 

ggplot( dat , aes(y= LD_APPROX , x = LD  )) + geom_point() 

 A: If you look at the theoretical approximation formula for LDi, you will see that LDi is a scalar of dimension 1x1 which is obtained as the product of three quantities: 


*

*The row vector of score residuals of dimension 1xp;

*The variance-covariance matrix of betahat of dimension pxp;

*The column vector of score residuals of dimension px1. 


I believe your R implementation of this formula should be:
t(SCORE) %% MOD$var %% SCORE
since R stores vectors as columns internally (despite the fact it prints them out in row format on the screen). There should be a star symbol included in between the double percentages above (but it gets removed for some reason when posting this).  Thus, you should not need to apply a diag() statement to the suggested R implementation, as it already produces a scalar (not a matrix). 
How good this approximation works in practice may depend on the peculiarities of your data. If you are not happy with it, just use the full expression for LDi (provided it is computed correctly).
