Semireplicated + REML in R [🇷 for BE/BA]

posted by Helmut Homepage – Vienna, Austria, 2020-07-10 21:03 (1771 d 06:45 ago) – Posting: # 21670
Views: 26,796

Hi ElMaestro,

❝ I am writing a little optimizer based on REML in R for use with semi-replicated designs. I am starting from the bottom and working upwards, since the various packages do not really help when R is replicated and T is not.


Wow! Nobody mastered that in R so far. I don’t want to dampen your enthusiasm but I think that is not possible at all.

❝ Does anyone have a dataset with well-established REML reference values for the variance components for the TRR/RRT/RTR design? A little more than the CVref of EMA dataset II would be appreciated :-)


Here you are (for DS II, of course). Results in Phoenix for the FDA’s covariance structure „Banded No-Diagonal Factor Analytic (f=2)” which is in SAS-lingo FA0(2). After 7 iterations (-2REML -48.2, AIC -26.2):
lambda(1,1) 0.19030501
lambda(1,2) 0.23981134
lambda(2,2) 0.074350441
s2wR        0.013246493
s2wT        0.0074570998
G Matrix*
0.03621600 0.04563730
0.04563730 0.06303747
As usual:
WARNING: Newton's algorithm converged with modified Hessian. Output is suspect.
Model may be over-specified. A simpler model could be tried.
s2wT is nonsense, since T is not replicated and the model is over-specified.
The lambdas are called in SAS-lingo FA(x,y). After 7 iterations (-2REML -30.7, AIC -20.7):
Cov Parm     Subject    Group    Estimate
FA(1,1)      Subj                  0.1903
FA(2,1)      Subj                  0.2398
FA(2,2)      Subj                  0.1072
Residual     Subj       trt R      0.01325
Residual     Subj       trt T      0.001497

    Estimated G Matrix
Row  trt    Col1      Col2
 1    R   0.03622   0.04564
 2    T   0.04564   0.06900
Also:
NOTE: Convergence criteria met but final hessian is not positive definite.
NOTE: Asymptotic variance matrix of covariance parameter estimates has been found to be
      singular and a generalized inverse was used. Covariance parameters with zero variance do
      not contribute to degrees of freedom computed by DDFM=SATTERTH.

No problems with convergence with f=1 (or FA0(1) in SAS). FA0(1) helps in general but not always.
lambda(1,1) 0.19030501
lambda(1,2) 0.23981134
s2wR        0.013246493
s2wT        0.012985088
Different s2wT but crap as well.

[image]You can get anything you like. The same data set throwing warnings in one software but not the other; same with convergence. In Certara’s forum a user posted a data set which failed to converge both in SAS and Phoenix, even with FA0(1). Then one is in deep shit.

Even if a data set converges without warnings, you might get different estimates (except the first one and s2wR) in SAS and Phoenix. Not surprising cause the optimizer walks around in parameter-space trying to find the maximum of a surface which is essentially flat. Increasing the number of iterations and/or decreasing the tolerance rarely helps.

I still wait for someone explaining why he/she is using such a wacky design. Beyond me.

—
Dif-tor heh smusma 🖖🏼 Довге життя Україна! [image]
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes

Complete thread:

UA Flag
Activity
 Admin contact
23,424 posts in 4,927 threads, 1,669 registered users;
182 visitors (0 registered, 182 guests [including 8 identified bots]).
Forum time: 03:49 CEST (Europe/Vienna)

By all means let’s be open-minded, but not so open-minded
that our brains drop out.    Richard Dawkins

The Bioequivalence and Bioavailability Forum is hosted by
BEBAC Ing. Helmut Schütz
HTML5