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

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

G Matrix*

As usual:

s2wT is nonsense, since T is not replicated and the model

The lambdas are called in SAS-lingo

Also:

No problems with convergence with f=1 (or

Different s2wT but crap as well.

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

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

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

» 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.

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

\(\begin{matrix}

s_{bT}^2 & \textrm{cov}(bT,bR)\\

\textrm{cov}(bT,bR) & s_{bR}^2

\end{matrix}\),

where \(s_{bT}^2\) and \(s_{bR}^2\) are the between subject variances of T and R, respectively and \(\textrm{cov}(bT,bR)\) is the between subject covariance for T and R given by \(s_{bT}\times s_{bR}\). The subject-by-formulation interaction variance is given by \(s_{bT}^2+s_{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_{bT}^2\) and \(s_{wT}^2\) cannot be unambiguously estimated, only their sum. Even then, PHX and SAS don’t agree (0.04367*vs.*0.03771).

—

Helmut Schütz

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

Science Quotes

*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 [R for BE/BA]
- 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

- "we"? Loosers? ElMaestro 2020-07-14 15:07

- Are we loosers? Helmut 2020-07-14 13:58

- took just 52 hrs to do it :-) ElMaestro 2020-07-13 14:18

- Medium rare. Helmut 2020-07-13 13:52

- 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

- Who can help? ElMaestro 2020-07-12 12:28

- Braveheart! Helmut 2020-07-11 15:38

- Avoid partial replicate designs, pleeeze! ElMaestro 2020-07-11 14:43

- Avoid partial replicate designs, pleeeze! Helmut 2020-07-11 11:57

- Semireplicated + REML in R ElMaestro 2020-07-10 19:24
- 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

- FDA RSABE is ISC d_labes 2020-07-15 18:13

- Desultory thoughts Helmut 2020-07-15 17:33
- Phoenix - which template? mittyri 2020-07-19 00:42
- FDA RSABE Project template_ v1.4.phxproj Helmut 2020-07-19 01:45

- It is the opposite way around for me ElMaestro 2020-07-15 16:25
- "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

- Speed improvement ElMaestro 2020-08-08 12:25

- Still can't make it work ElMaestro 2020-08-07 18:23

- Still can't make it work PharmCat 2020-08-07 18:14

- Still can't make it work ElMaestro 2020-08-07 16:44

- Still can't make it work PharmCat 2020-08-07 15:42

- Still can't make it work ElMaestro 2020-08-07 13:08

- Open issues PharmCat 2020-08-07 11:29

- Open issues ElMaestro 2020-08-07 07:49

- Open issues PharmCat 2020-08-07 00:02

- Open issues ElMaestro 2020-08-06 21:11

- Not understood PharmCat 2020-08-05 16:37

- Not understood ElMaestro 2020-08-05 08:13

- Not understood PharmCat 2020-08-05 01:41

- Not understood ElMaestro 2020-08-03 22:55

- And by the way.... PharmCat 2020-08-03 14:24

- And by the way.... ElMaestro 2020-07-24 12:52

- Semireplicated + REML in RHelmut 2020-07-10 19:03