06-05-2014 05:20 AM
I simulated a dataset with 500 profiles of a biomarker measured every 3 weeks during 2 years. My structural model is physiological and so a little complex. It can be compared to a biexponentiel function with a time of rupture D. I have 4 fixed effects, 4 random effects (3 are log-transformed and 1 logit-transformed) and a constant residual error on log(Y+1) to estimate. My code is the following:
proc nlmixed data=dataset;
parms mu_A=0.05 mu_B=80 mu_C=0.3 mu_D=140
omega_A=0.1 omega_B=0.6 omega_C=1.5 omega_D=0.6
if (TIME<D) then Y= 0.23*B/(A*(1-C)-0.046+0.23)*exp((A*(1-C)-0.046)*TIME)+(B-0.23*B/(A*(1-C)-0.046+0.23))*exp(-0.23*TIME) ;
else Y=0.23*B*exp((A*(1-C)-0.046)*D) /(A-0.046+0.23)*exp((A-0.046)*(TIME-D))+(YD-0.23*B*exp((A*(1-C)-0.046)*D)/(A-0.046+0.23))*exp(-0.23*(TIME-D)) ;
model logY ~ normal(lY, sigma**2);
random eta_A eta_B eta_C eta_D ~ normal([0,0,0,0], [omega_A**2,
0, 0, omega_C**2,
0, 0, 0, omega_D**2]) subject=ID;
I would like to use the adaptative Gauss quadrature or Laplace approximation. This code works for the methods FIRO or for non-adaptative Gauss quadrature (but NOAD needs a lot of QPOINTS and takes a long time). But for the adaptative Gauss quadrature, I have the error message: No valid parameter points were found and the initial negative log-likelihood is 1.15792E77. I don’t understand this error since my dataset is simulated given this model and I initialize parameters with the good values. I tried simpler models, notably with only one or two random effects or with smaller omega. Some work, if omegas are enough small, or if only A and B have random effects.
Is my model too complex for SAS? Is there a problem with the rupture time D? Have I too much variability?
My goal is then to write my own likelihood function and to fit real data, so I could not control their variability.
Thanks in advance
06-06-2014 09:20 AM
Adaptive quadrature gets problematic with multiple random effects and a structured covariance AND complex entry into the model. That is why the NOAD works, although slowly. The variance for C looks like the real source of the blow-up, but I don't see a good re-parameterization from the logistic. Could you "unstructure" the covariance matrix (this means adding six more parameters to optimize), but constrain them somehow.
Check for typos--I see D shows as YD at least once in the else clause.
Consider rescaling time. The exponentiating can lead to problems. Another approach would be to model log(time).
08-03-2016 05:45 AM
What I have tried and it helped is that trying to use "good starting values".
I also have simulation data. My model had one RE. And I used the starting values that were quite close to the "true value" (value used in the simulation) but it produced the error "No valid parameter points were found".
1. Remove the RE, fit the model with initial values qutie close to the true value.
The model converged well. I got the paramter estimates.
2. I used paramter estimates from step 1 as starting values for model with RE, now the initial value for the variance of RE is the value close to one that I used in simulation step.
The model run, but now with warning of "Hessian matrix is full rank but has at least one negative eigenvalue".
I try to vary the initial value of the variance of the RE and finally I got the results without any warnings.