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

posted by Helmut Homepage – Vienna, Austria, 2020-07-10 19:03 (82 d 10:45 ago) – Posting: # 21670
Views: 10,276

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

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 🖖
Helmut Schütz

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

Complete thread:

 Admin contact
21,090 posts in 4,398 threads, 1,469 registered users;
online 13 (0 registered, 13 guests [including 5 identified bots]).
Forum time: Thursday 05:49 CEST (Europe/Vienna)

In these days, a man who says a thing cannot be done
is quite apt to be interrupted by some idiot doing it.    Elbert Green Hubbard

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