Proc IML -> Verifying -2 Log-Likelihood Issue; AR(1) Matrix, CS Matrix

[mattll218]

Hi all!

Forgive me in advance, I am in no way a SAS professional and am trying to learn. I'm working on version 9.4 for reference.

I am working on a code to compute the estimated log-likelihood for a model given parameters, particularly to verify the maximized log-likelihood function value. The covariance is being analyzed with an AR (1) matrix and a CS matrix. The computed log-likelihood using proc mixed for the AR (1) matrix is -256.7323 and CS matrix is -204.7334, but the following code throws log-likelihoods that are extremely off. Part a.2-3 is for the AR (1) matrix and Part A.4 is for the CS Matrix. There is no missing data in the dataset. If anyone could point me to a resource or has a suggestion, please let me know!

/*Part A.2-3*/
proc iml;
   start log_normal_density(y, mu, sigma, rho);
       n = nrow(y);
       R = j(n, n, .);
       do i = 1 to n;
           R[i,i] = 1;
           do j = i+1 to n;
               R[i,j] = rho##abs(i-j);
               R[j,i] = R[i,j];
           end;
       end;
       sigmaMat = sigma##2 * R;
       eigenvalues = eigval(sigmaMat);
       logDetSigma = sum(log(eigenvalues));
       if logDetSigma <= 0 then logDetSigma = 1e-8;
       invSigmaMat = inv(sigmaMat);
       const = -0.5 * n * log(2 * constant("pi")) - 0.5 * logDetSigma;
       diff = y - mu;
       quad_form = sum(diff # invSigmaMat * diff);
       log_pdf = const - 0.5 * quad_form;
       return(log_pdf);
   finish;
   start mylik(y, a, b, sigma, rho);
       n = nrow(y);
       mu = (a + b * (0:n-1))`;
       lik = log_normal_density(y, mu, sigma, rho);
       return(lik);
   finish;
   use rats;
   read all var {wei} into Y;
   close rats;
   a = 4.0237;
   b = 0.2458;
   sigma = sqrt(0.01764);
   rho = 0.7643;
   logL = mylik(Y, a, b, sigma, rho);
   neg2LogL = (-2) * logL;
   print logL neg2LogL;
quit;
/*Part A.4*/
proc iml;
start log_normal_density_cs(y, mu, sigma, rho);
    n = nrow(y);

    /* Correct CS covariance matrix */
    R = j(n, n, rho);        
    do i = 1 to n;
        R[i,i] = 1;    
    end;

    sigmaMat = sigma * ((1 - rho) * I(n) + rho * J(n, n, 1)); 

    detR = det(sigmaMat);    
    invR = inv(sigmaMat);     

    const = -0.5 * n * log(2 * constant("pi")) - 0.5 * log(detR);
    log_pdf = const - 0.5 * t(y - mu) * invR * (y - mu);

    return(log_pdf);
finish;

start mylik_cs(y, a, b, sigma, rho);
    n = nrow(y);
    mu = a + b * T(0:n-1);  /* Transpose to match y's dimension */
    lik = log_normal_density_cs(y, mu, sigma, rho);
    return(lik);
finish;
use rats;
read all var {wei} into Y;
close rats;
a = 4.0615;
b = 0.2431;
sigma = 0.009249;
rho = 0.007725;
logL_cs = mylik_cs(Y, a, b, sigma, rho);
neg2LogL_cs = (-2) * logL_cs;

print logL_cs neg2LogL_cs;
quit;

[Rick_SAS]

I'll try to look at your code more closely later today, but here are some initial suggestions:

1. You can get the AR(1) correlation matrix by using the TOEPLITZ function. See the AR1COEFF function defined in Two simple ways to construct a log-likelihood function in SAS - The DO Loop

2. If you want to compute the log-determinant of a matrix, use the LOGABSDET function. See Compute the log-determinant of an arbitrary matrix - The DO Loop

3. To compute the log-density of common probability distributions, use the LOGPDF function. Again, see Two simple ways to construct a log-likelihood function in SAS - The DO Loop

Can you post the RATS data?

[Rick_SAS]

What does your PROC MIXED code look like? Can you upload the data?