## Still can't make it work [R for BE/BA]

Hi PharmCat.

» If you calculate β with initial V - it will be wrong.

Correct and that is not what I am doing. For every iteration through the optimizer I am deriving a new beta vector; this is beta given V. The loglikelihood is calculated on basis of V and beta, where beta is derived from V via X. Therefore LL can be said to be a function of V (or of the variance components in it).

There are two alternative ways I can get that beta vector for any given V:

1. I can do it the simple way with the very big X. It works well.

2. I can perhaps (=I should be able to?) do it by partition X into small bits, Xi, like in Gurka equation 3, or Lindstrom&Bates. But I cannot get the same beta vector for a given set of variance components in V.

The -or a- reason is that if I reduce X to Xi (individual subjects) and do all the summation then the M1 and written above is not invertible. This, in turn, appears to be because the individual Xi's are not of the same rank as X. If I lump Xi's together like the first 6 subjects, then the next 6 subjects etc, and so forth and use Gurka's equation 3, then I can generate a beta estimate (from 4 chunks of 6 subjects),

I am not even at a final iteration, I am just at any iteration where I am comparing beta via the two methods for a set of starting values, or for that matter, if I plug in the values that come out as optimal values from pt 1 above after everything is coverged.

» It will be estimate for current V, then you calk REML end opimize θ... And then again... recalk β, ets... You can get initial β as β = inv(qr(X).R)*qr(X).Q'*y

You do not need this since you do not need starting values for beta, cf the LL above. You only need starter values for theta (variance components in V). Every pass through the REML LL equation generates a beta estimate from the variance compnents, given X. Once the algo has converge d the estimate of beta is final.

All this imaterial to my challenge at hand.

» But if you have wrong β at final iteration - it is bad case.

This is not about the "final iterattion". It is about any iteration. For any iteraiton, the two ways to derive beta results in different values. Therefore something s wrong.

» You can find same equation in LINDSTROM&BATES - 1.2 (where X not a design matix, but vector of X

» If you provide me data from 6 subjects I can make calculations in Julia.

Here all subjects are (y is logarithmised):

Subj Seq Trt Per y

1 RTR R 1 8.307361

1 RTR T 2 8.286622

1 RTR R 3 8.229191

2 RTR R 1 8.001757

2 RTR T 2 7.774351

2 RTR R 3 7.939016

3 RRT R 1 8.150295

3 RRT R 2 8.113786

3 RRT T 3 8.301224

4 TRR T 1 8.319961

4 TRR R 2 8.068152

4 TRR R 3 8.243703

5 TRR T 1 8.469640

5 TRR R 2 8.421255

5 TRR R 3 8.278936

6 RRT R 1 8.023159

6 RRT R 2 8.015393

6 RRT T 3 7.791358

7 RTR R 1 7.836054

7 RTR T 2 8.030084

7 RTR R 3 7.993823

8 TRR T 1 7.698483

8 TRR R 2 7.621391

8 TRR R 3 7.609862

9 RTR R 1 8.444106

9 RTR T 2 8.333174

9 RTR R 3 8.131531

10 TRR T 1 7.708949

10 TRR R 2 7.766586

10 TRR R 3 7.705803

11 RRT R 1 7.530373

11 RRT R 2 7.701833

11 RRT T 3 7.780888

12 RRT R 1 7.731229

12 RRT R 2 8.061613

12 RRT T 3 8.275682

13 RRT R 1 7.878686

13 RRT R 2 7.795811

13 RRT T 3 7.961789

14 RTR R 1 8.016582

14 RTR T 2 7.807835

14 RTR R 3 7.996452

15 RTR R 1 7.720639

15 RTR T 2 7.598299

15 RTR R 3 7.910003

16 TRR T 1 7.992809

16 TRR R 2 8.143808

16 TRR R 3 8.114504

17 TRR T 1 7.781890

17 TRR R 2 7.885856

17 TRR R 3 7.683404

18 RRT R 1 7.910224

18 RRT R 2 7.939373

18 RRT T 3 8.054078

19 RTR R 1 7.791027

19 RTR T 2 7.919283

19 RTR R 3 7.825645

20 RRT R 1 7.886983

20 RRT R 2 7.982689

20 RRT T 3 8.018691

21 RRT R 1 7.961928

21 RRT R 2 7.888485

21 RRT T 3 8.029107

22 TRR T 1 7.989221

22 TRR R 2 7.981665

22 TRR R 3 7.956967

23 TRR T 1 8.056680

23 TRR R 2 8.066396

23 TRR R 3 8.174308

24 RTR R 1 7.536257

24 RTR T 2 7.500419

24 RTR R 3 7.912350

» If you calculate β with initial V - it will be wrong.

Correct and that is not what I am doing. For every iteration through the optimizer I am deriving a new beta vector; this is beta given V. The loglikelihood is calculated on basis of V and beta, where beta is derived from V via X. Therefore LL can be said to be a function of V (or of the variance components in it).

There are two alternative ways I can get that beta vector for any given V:

1. I can do it the simple way with the very big X. It works well.

2. I can perhaps (=I should be able to?) do it by partition X into small bits, Xi, like in Gurka equation 3, or Lindstrom&Bates. But I cannot get the same beta vector for a given set of variance components in V.

The -or a- reason is that if I reduce X to Xi (individual subjects) and do all the summation then the M1 and written above is not invertible. This, in turn, appears to be because the individual Xi's are not of the same rank as X. If I lump Xi's together like the first 6 subjects, then the next 6 subjects etc, and so forth and use Gurka's equation 3, then I can generate a beta estimate (from 4 chunks of 6 subjects),

*but it is still not the same as in pt 1*, for any given set of components in V.I am not even at a final iteration, I am just at any iteration where I am comparing beta via the two methods for a set of starting values, or for that matter, if I plug in the values that come out as optimal values from pt 1 above after everything is coverged.

» It will be estimate for current V, then you calk REML end opimize θ... And then again... recalk β, ets... You can get initial β as β = inv(qr(X).R)*qr(X).Q'*y

You do not need this since you do not need starting values for beta, cf the LL above. You only need starter values for theta (variance components in V). Every pass through the REML LL equation generates a beta estimate from the variance compnents, given X. Once the algo has converge d the estimate of beta is final.

All this imaterial to my challenge at hand.

» But if you have wrong β at final iteration - it is bad case.

This is not about the "final iterattion". It is about any iteration. For any iteraiton, the two ways to derive beta results in different values. Therefore something s wrong.

» You can find same equation in LINDSTROM&BATES - 1.2 (where X not a design matix, but vector of X

_{i}, IMHO)» If you provide me data from 6 subjects I can make calculations in Julia.

Here all subjects are (y is logarithmised):

Subj Seq Trt Per y

1 RTR R 1 8.307361

1 RTR T 2 8.286622

1 RTR R 3 8.229191

2 RTR R 1 8.001757

2 RTR T 2 7.774351

2 RTR R 3 7.939016

3 RRT R 1 8.150295

3 RRT R 2 8.113786

3 RRT T 3 8.301224

4 TRR T 1 8.319961

4 TRR R 2 8.068152

4 TRR R 3 8.243703

5 TRR T 1 8.469640

5 TRR R 2 8.421255

5 TRR R 3 8.278936

6 RRT R 1 8.023159

6 RRT R 2 8.015393

6 RRT T 3 7.791358

7 RTR R 1 7.836054

7 RTR T 2 8.030084

7 RTR R 3 7.993823

8 TRR T 1 7.698483

8 TRR R 2 7.621391

8 TRR R 3 7.609862

9 RTR R 1 8.444106

9 RTR T 2 8.333174

9 RTR R 3 8.131531

10 TRR T 1 7.708949

10 TRR R 2 7.766586

10 TRR R 3 7.705803

11 RRT R 1 7.530373

11 RRT R 2 7.701833

11 RRT T 3 7.780888

12 RRT R 1 7.731229

12 RRT R 2 8.061613

12 RRT T 3 8.275682

13 RRT R 1 7.878686

13 RRT R 2 7.795811

13 RRT T 3 7.961789

14 RTR R 1 8.016582

14 RTR T 2 7.807835

14 RTR R 3 7.996452

15 RTR R 1 7.720639

15 RTR T 2 7.598299

15 RTR R 3 7.910003

16 TRR T 1 7.992809

16 TRR R 2 8.143808

16 TRR R 3 8.114504

17 TRR T 1 7.781890

17 TRR R 2 7.885856

17 TRR R 3 7.683404

18 RRT R 1 7.910224

18 RRT R 2 7.939373

18 RRT T 3 8.054078

19 RTR R 1 7.791027

19 RTR T 2 7.919283

19 RTR R 3 7.825645

20 RRT R 1 7.886983

20 RRT R 2 7.982689

20 RRT T 3 8.018691

21 RRT R 1 7.961928

21 RRT R 2 7.888485

21 RRT T 3 8.029107

22 TRR T 1 7.989221

22 TRR R 2 7.981665

22 TRR R 3 7.956967

23 TRR T 1 8.056680

23 TRR R 2 8.066396

23 TRR R 3 8.174308

24 RTR R 1 7.536257

24 RTR T 2 7.500419

24 RTR R 3 7.912350

—

I could be wrong, but...

Best regards,

ElMaestro

R's base package has 274 reserved words and operators, along with 1761 functions. I can use 18 of them (about 14 of them properly). I believe this makes me the Donald Trump of programming.

I could be wrong, but...

Best regards,

ElMaestro

R's base package has 274 reserved words and operators, along with 1761 functions. I can use 18 of them (about 14 of them properly). I believe this makes me the Donald Trump of programming.

### Complete thread:

- Semireplicated + REML in R ElMaestro 2020-07-10 13:19 [R for BE/BA]
- Semireplicated + REML in R Helmut 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 workElMaestro 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 workElMaestro 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 R Helmut 2020-07-10 19:03