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

posted by ElMaestro  – Belgium?, 2020-08-07 16:44 (48 d 12:04 ago) – Posting: # 21834
Views: 5,931

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), 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 Xi, 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.

Complete thread:

Activity
 Admin contact
21,074 posts in 4,394 threads, 1,468 registered users;
online 2 (0 registered, 2 guests [including 2 identified bots]).
Forum time: Friday 04:49 CEST (Europe/Vienna)

If you think it’s simple,
then you have misunderstood the problem.    Bjarne Stroustrup

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