# Optimize likelihood function to get lambda for Box-Cox transform of two variables

I'm trying to implement the Transformed Normal method (TN) that was described in Faraggi, D., & Reiser, B. (2002). Estimation of the area under the ROC curve. Statistics in medicine, 21(20), 3093–3106.

The authors transform the measurements of a diseased and a non-diseased sample to normality using a common $\hat{\lambda}$ in the Box-Cox-transformation and then proceed to calculate a cutpoint using normal theory.

The $\hat{\lambda}$ is obtained according to Zou, K. H., Tempany, C. M., Fielding, J. R., & Silverman, S. G. (1998). Original smooth receiver operating characteristic curve estimation from continuous data: Statistical methods for analyzing the predictive value of spiral CT of ureteral stones. Academic Radiology, 5(10), 680–687. https://doi.org/10.1016/S1076-6332(98)80562-X

Here's the relevant part of the latter paper: So a profile log likelihood function is constructed and optimized to find $\hat{\lambda}$. $x_{0,i}$ and $y_{0,j}$ are the original, unaltered measurements and $x'$ and $y'$ are the Box-Cox-transformed ones, if I understand the paper correctly.

I am now trying to implement this likelihood function in R. Apparently, the Box-Cox-transformed variables are inserted into the likelihood function, so my attempt looks like this:

boxcox <- function(x, lambda) {
if (lambda == 0) {
return(log(x))
} else {
((x**lambda) - 1) / lambda
}
}

profile_loglik <- function(lambda, x, y) {
m <- length(x)
n <- length(y)
x_t <- boxcox(x, lambda)
y_t <- boxcox(y, lambda)
return(-m * log(sd(x_t)) - n * log(sd(y_t)) + (lambda - 1) * (sum(log(x_t)) + sum(log(y_t))))
}

calc_lambda_hat <- function(x, y) {
optim_func <- function(lambda) -(profile_loglik(lambda = lambda, x = x, y = y))
optim(0.5, optim_func, method = "Brent", lower = 0, upper = 100)
}

# Example data
x <- rlnorm(200)
y <- rlnorm(100, 0.4, 0.4)
calc_lambda_hat(x, y)

# Result:
$par  100$value
 NaN

$counts function gradient NA NA$convergence
 0