Bioequivalence and Bioavailability Forum

In matrix notation will write

y=Xb+e

y is a vector of observations.
X is the design matrix.
b is the effects vector.
e is the error or residual vector.

Let's somehow try to get the dataset above transmogriffed into some kind of matrix notation.

I will not go through matrix notation or mutiplication here. Google it, or use Wikipedia. Along the same lines, I will not go through anything that leads to the least squares solution of the normal linear model, but I will mention without proof that:

1. If the error is normal with mean zero and some variance then the most likely set of constants in b is those for which the sum of squared residuals ete is minimal.
   That solution is analytically
   b=(XtX))-1Xty
   
   
2. This solution implies inversion of the matrix XtX. If this matrix cannot be inverted then we have screwed up.
   
   
3. For all practical purposes, if XtX is invertible, then we have a solution. This is the case when the rank of X equals the number of columns in X.
   
   
4. If any columns of X add up to give another column (or sum of columns) in X then the rank is not equal to the number of columns, and we cannot get a solution.

We can write the equations like this (focus for now only on sequence and observations):

17.11=1*intercept + 1*Seq_TR + 0*Seq_RT + …
17.76=1*intercept + 1*Seq_TR + 0*Seq_RT + …
…
…
…
22.87=1*intercept + 0*Seq_TR + 1*Seq_RT + …
20.80=1*intercept + 0*Seq_TR + 1*Seq_RT + …


Let's transmogriph this into a model matrix. In X we have a column for each constant we wish to determine. Above we just worked with three constants (intercept, and two sequence constants).
In the rows of the intercept put a one if the observations is associated with an intercept. All observations are so this is a column of pure 1's.
In the rows of Sequence TR put a 1 when the observation comes from a subject in sequence TR otherwise a 0. And so forth for Sequence RT.

We can write the model matrix (shown right in reddish color) like:



Note the bold constants just above: They correspond 1:1 to how they appear in X.

But there's trouble ahead.

Look at column 2 and 3 of X (for sequence TR and RT, respectively). If you add them, then you get a column of pure 1's which is the same as column 1 for the intercept.

In other words, X writen just above does not satisfy the rules for finding a solution, cf. rules above. Take away e.g. the column for SeqRT and then you can see that addition of one column to some other column(s) does not result in some other column(s). When we've taken away the SeqRT column, then the rank is therefore two, and we have two columns, and therefore X^tX is invertible, and therefore we have a solution.
In other words, We have two sequences but we can only fit one constant for sequence (a constant for either SeqTR or SeqRT, assuming we insist on having an intercept). In other words: We just lost a degree of freedom. And we are not getting it back.

So, the unfinished X now looks like this:




Now let's add two columns for period just using the same principles (X to the right):



And again, this X is sick, ill and syphilitic because the two period columns add up to the intercept. We need to get rid of one of these period columns and then all is good again.



Now XtX is invertible and we're in the clear.

We do the same for treatment (can't have two columns for treatment since they add up to the intercept) and get:

Now all that is left is to add columns for subject:



(first two subject columns shown in X)

If we have N subjects then we'd initially try to add N columns to X, but it turns out that we need to eliminate 2 of the subject columns in order to keep X^tX invertible (corresponding to a full rank X) so that we can find the maximum likelihood solution.

The status is

• We wanted a model with an intercept and four fixed factors called Sequence, Period, Treatment, Subject.
  
• Intercept has 1 df.
  
• We lost one DF for sequence. Sequence has 1 df.
  
• We lost one DF for period. Period has 1 df.
  
• We lost one DF for treatment. Treatment has 1 df.
  
• We lost two DFs for subject. Subject has N-2 df's.

The vector b will thus have one row for each column in X, in total there will be 1+1+1+1+(N-2) rows in b.

Nifty stuff: If you know the subject, the period, the treatment for any observation then you know tha subject's sequence. And so forth. This is why the way I wrote the model in the beginning of this thread can be smartified like C&L did it. ~ A loss of a df means something can be figured out if we just have access to auxiliary info.


But come on...
We tried to get two columns for Sequence but we had to remove one due to the intercept. Same for Period and Treatment. For Subject it got even worse. We can just fit the model without intercept, right? Then all trouble is gone!

Sort of. Then we can add two columns for Sequence (if this is the first term we deal with), now the two sequence columns add up to pure 1's, so we cannot add two period columns, and we cannot add two treatment columns, but just one of each again.

All is good so far.

But here's the bad: The dark forces are upon us. The lurk on every corner and we cannot escape the consequences of their hideous actions. They talk about... ... oh, forgive the typos, for my hands are shaking the very moment I think of it... they talk about "Subject nested in Sequence" and feel smug about it because it sounds fancy.


AAAAAAAAAAaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaarggggh. Gimme some garlic and some holy water. Make them go away, send them to work in the salt mines for 50 years, anything. ANYTHING!


When we talk about "Subject nested in Sequence" or "Subject in Sequence" or "Subject(Sequence)", which is all the same, then it means we try to introduce columns in X corresponding to any combo of Subject (N levels) AND Sequence (2 levels).
Let's just look at a few columns in X for subject 1 in sequence TR and subject 1 in sequence RT:


See? This is the ugly. There was only one subject called "Subject 1", she was assigned to sequence TR, so it makes absolutely no sense to speak of subject 1 in sequence RT; the corresponding column will be one of pure zero's which does not convey any info of any sort (if 0*a = 0 then what can you tell me about a?).

Does it means that people who talk about "Subject nested in Sequence" are full of sh!t?
If you ask me, then by and large yes. In some cases it makes a little sense, though.
Imagine we do not call our subjects 1,2,3 ... N like in the examples above.

In stead we assign number 1,2,3... to each subject in sequence TR.
And we assign number 1,2,3... to each subject in sequence RT.
And we make sure to keep those numbers apart. Talking about "subject 1" is now no longer sufficent for identification, so to identify the subject we need to know which sequence he/she was assigned to.

If you number your subjects this way and you construct your X following the rules above then it will make a little sense to do it, and you will arrive at EXACTLY the same columns in X as if you do it with unquely numbered subjects.

In my years of working with clinical documentation I have never ever seen a CRO or clinic use anything but uniquely numbered subject identifiers. Actually these are always the randomsation codes and these are always unique, right? This is why it really doesn't make practical sense to speak of "Subject nested in Sequence". It isn't wrong per se but it involves a lot ridiculous and meaningless matrix algebra which doesn't give any info at all.
So stay good: if you work with unique randomisation numbers assigned to the subjects (and you do, I know you do!) then just forget about the nesting stuff.

So we had a dataset, we got X constructed so that it was full rank and we didn't waste any time on meaningless babble about subjects nested in sequence. The number of columns of X is the rank, which is the number of rows in b, which is the number of constants we have freedom to determine. We're doing this by minimising ete. The story ends here, submit, get approval, celebrate.

Almost. I often do not fit with the intercept, and then I specify Treatment as the first factor. This means we do not have an interceept column in X, and the first two columns will be the two treatment columns. This is very handy, because then the first two estimates or constants in b are the treatment effects: One for Test and one for Reference.
So we just extract these to estimates from the b-vector and use them directly for the construction of the confidence interval. Note that the residual sum of squares is not affected by this. It doesn't matter which factors you specify first or whether you want an intercept. If you have one way or another Period, Subject, Sequence, Treatment in your model, then regardless of the permutation of the order of them the mean squared error is unaffected.

The question arises if the English expression and comprehension is more difficult than the MATHEMATICS. Now for the 12 year old.

My understanding is the mixed model (in my case from Phoenix WinNonlin) predicts the value of the PK parameter by a regression model in the terms that you describe. You then can compare with the observed value of the parameter.

Now I think this is correct, but by the way why do we call it a mixed model? I am trying to shape my understanding of this activity better

ANGUS