Bioequivalence and Bioavailability Forum • Semireplicated + REML in R

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

posted by Helmut – Vienna, Austria, 2020-07-10 21:03 (1378 d 21:27 ago) – Posting: # 21670
Views: 18,277

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.

The G (variance-covariance) matrix is obtained by
\(\small{\begin{matrix}
s_\textrm{bT}^2 & \textrm{cov(bT,bR)}\\
\textrm{cov(bT,bR)} & s_\textrm{bR}^2
\end{matrix}}\),
where \(\small{s_\textrm{bT}^2}\) and \(\small{s_\textrm{bR}^2}\) are the between subject variances of T and R, respectively and \(\small{\textrm{cov(bT,bR)}}\) is the between subject covariance for T and R given by \(\small{s_\textrm{bT}\times s_\textrm{bR}}\). The subject-by-formulation interaction variance is given by \(\small{s_\textrm{bT}^2+s_\textrm{bR}^2-2\times\textrm{cov(bT,bR)}}.\)
Whereas one might get a negative lambda(2,2) in PHX, in such a case SAS forces FA(2,2) to zero. With both approaches one might end up with a negative subject-by-formulation interaction variance…
The same holds for the within subject variance of T (can be negative in PHX but is forced to zero in SAS).
However, in partial replicate designs \(s_\textrm{bT}^2\) and \(s_\textrm{wT}^2\) cannot be unambiguously estimated, only their sum. Even then, PHX and SAS don’t agree (0.04367 vs. 0.03771).

—
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:

Semireplicated + REML in R ElMaestro 2020-07-10 13:19
- Semireplicated + REML in RHelmut 2020-07-10 19:03
  - Semireplicated + REML in R ElMaestro 2020-07-10 19:24
    - Avoid partial replicate designs, pleeeze! Helmut 2020-07-11 11:57
      - Avoid partial replicate designs, pleeeze! ElMaestro 2020-07-11 14:43
        
        Braveheart! Helmut 2020-07-11 15:38
        
        Who can help? ElMaestro 2020-07-12 12:28
        
        Who can help? ElMaestro 2020-07-12 13:26
        Update II ElMaestro 2020-07-12 21:36
        
        Update III ElMaestro 2020-07-12 21:46
        
        Final update today ElMaestro 2020-07-12 22:27
        
        Medium rare. Helmut 2020-07-13 13:52
        
        took just 52 hrs to do it :-) ElMaestro 2020-07-13 14:18
        
        Will take much more hours still… Helmut 2020-07-13 15:34
        
        Negative determinant ElMaestro 2020-07-14 03:22
        
        Are we loosers? Helmut 2020-07-14 13:58
        
        "we"? Loosers? ElMaestro 2020-07-14 15:07
        
        Misunderstanding? Helmut 2020-07-14 15:32
        
        Braveheart! ElMaestro 2020-07-13 10:13
        
        Braveheart! Helmut 2020-07-13 14:16
        Braveheart! PharmCat 2020-07-15 14:19
        
        Braveheart! ElMaestro 2020-08-02 17:39
- We were all blind (except Detlew) Helmut 2020-07-15 14:27
  - It is the opposite way around for me ElMaestro 2020-07-15 16:25
    - Desultory thoughts Helmut 2020-07-15 17:33
      - FDA RSABE is ISC d_labes 2020-07-15 18:13
        
        FDA RSABE is ISC Helmut 2020-07-16 11:11
      - Desultory thoughts ElMaestro 2020-07-15 23:06
        
        Desultory thoughts Helmut 2020-07-16 10:59
  - Phoenix - which template? mittyri 2020-07-19 00:42
    - FDA RSABE Project template_ v1.4.phxproj Helmut 2020-07-19 01:45
- "By popular demand": likelihood ElMaestro 2020-07-24 10:07
  - And by the way.... ElMaestro 2020-07-24 12:52
    - And by the way.... PharmCat 2020-08-03 14:24
      - Not understood ElMaestro 2020-08-03 22:55
        
        Not understood PharmCat 2020-08-05 01:41
        
        Not understood ElMaestro 2020-08-05 08:13
        
        Not understood PharmCat 2020-08-05 16:37
        
        Open issues ElMaestro 2020-08-06 21:11
        
        Open issues PharmCat 2020-08-07 00:02
        
        Open issues ElMaestro 2020-08-07 07:49
        
        Open issues PharmCat 2020-08-07 11:29
        
        Still can't make it work ElMaestro 2020-08-07 13:08
        
        Still can't make it work PharmCat 2020-08-07 15:42
        
        Still can't make it work ElMaestro 2020-08-07 16:44
        
        Still can't make it work PharmCat 2020-08-07 18:14
        
        Still can't make it work ElMaestro 2020-08-07 18:23
        
        And now it works ElMaestro 2020-08-07 21:31
        Still can't make it work PharmCat 2020-08-08 01:08
        
        Speed improvement ElMaestro 2020-08-08 12:25
        
        Speed improvement PharmCat 2020-08-08 17:27
        Speed improvement ElMaestro 2020-08-08 18:10
        
        Speed improvement PharmCat 2020-08-09 18:22
        Some tests... PharmCat 2020-08-10 11:48