## 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 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
$$\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