Bioequivalence and Bioavailability Forum • Parallel design [unequal variances] in Phoenix/WinNonlin

Helmut
★★★

Vienna, Austria,
2013-07-04 22:22
(4309 d 18:38 ago)

Posting: # 10938
Views: 14,435

Parallel design [unequal variances] in Phoenix/WinNonlin [Software]

Dear all,

some veterans might remember this thread from 2007 (!) where I showed that WinNonlin assumes equal variances in parallel designs.
This week Linda Hughes of Pharsight posted a coding at the Extranet which works both in Phoenix and ‘classical’ WinNonlin (I tested PHX 6.3 and WNL 5.3). Her solution follows a hint given at SAS’ knowledge base. Of course, our SAS-guru Detlew knew that already. ;-)

Here is my setup in [image]

(AFAIK the same one is coded in bear). First the original balanced data set with similar variances (T 0.1796, R 0.1292) and then a modified imbalanced data set with unequal variances (T 0.5639, R 0.1292).

disply <- function(table) {

    cat("\nTest/Reference with 95% Confidence Limits (90% CI)", fill=T)

    round(table, 2) }



# original data

PK <- ("

  AUC, treatment

  44.1, R

  33.6, R

  45.5, R

  35.3, R

  26.0, R

  38.2, R

  25.6, R

  58.0, R

  47.2, R

  17.5, R

  51.7, R

  24.5, R

  19.5, T

  67.2, T

  25.7, T

  33.6, T

  25.1, T

  44.1, T

  16.5, T

  47.3, T

  22.6, T

  36.3, T

  29.4, T

  18.3, T")

txtcon <- textConnection(PK)

PKdata <- read.table(txtcon, header=T, sep=",", strip.white=T)

close(txtcon)

PKdata$logAUC <- log(PKdata$AUC)

attach(PKdata)

relevel(treatment, ref = "R")

t1 <- t.test(logAUC ~ treatment, data=PKdata, var.equal=T, conf.level=0.90)

tbl1 <- matrix(

    c(as.numeric(100*exp(diff(t1$estimate))), sort(as.numeric(100*exp(-t1$conf.int)))),

    byrow=T, nrow=1)

dimnames(tbl1) <- list("t-test (equ. var.) ", c("Point Estimate", "Lower CL", "Upper CL"))

t2 <- t.test(logAUC ~ treatment, data=PKdata, var.equal=F, conf.level=0.90)

tbl2 <- matrix(

    c(as.numeric(100*exp(diff(t2$estimate))), sort(as.numeric(100*exp(-t2$conf.int)))),

    byrow=T, nrow=1)

dimnames(tbl2) <- list("Welch-Satterthwaite", c("Point Estimate", "Lower CL", "Upper CL"))

disply(tbl1)

disply(tbl2)



# modified data

PK <- ("

  AUC, treatment

  44.1, R

  33.6, R

  45.5, R

  35.3, R

  26.0, R

  38.2, R

  25.6, R

  58.0, R

  47.2, R

  17.5, R

  51.7, R

  24.5, R

  58.5, T

  201.6, T

  77.1, T

  33.6, T

  25.1, T

  44.1, T

  16.5, T

  47.3, T

  22.6, T")

txtcon <- textConnection(PK)

PKdata <- read.table(txtcon, header=T, sep=",", strip.white=T)

close(txtcon)

PKdata$logAUC <- log(PKdata$AUC)

attach(PKdata)

relevel(treatment, ref = "R")

t3 <- t.test(logAUC ~ treatment, data=PKdata, var.equal=T, conf.level=0.90)

tbl3 <- matrix(

    c(as.numeric(100*exp(diff(t3$estimate))), sort(as.numeric(100*exp(-t3$conf.int)))),

    byrow=T, nrow=1)

dimnames(tbl3) <- list("t-test (equ. var.) ", c("Point Estimate", "Lower CL", "Upper CL"))

t4 <- t.test(logAUC ~ treatment, data=PKdata, var.equal=F, conf.level=0.90)

tbl4 <- matrix(

    c(as.numeric(100*exp(diff(t4$estimate))), sort(as.numeric(100*exp(-t4$conf.int)))),

    byrow=T, nrow=1)

dimnames(tbl4) <- list("Welch-Satterthwaite", c("Point Estimate", "Lower CL", "Upper CL"))

disply(tbl3)

disply(tbl4)

In Phoenix you need a Subject variable and a dummy, which contains 1 in all cells. For simplicity I named it Period.
In the default coding you have only one Fixed Effect, namely Treatment. Now for the trick:
In the Setup map additionally Subject and Period to Classification.
Fixed: Treatment

Variance Structure > Repeated:

Repeated Specification: Period

()

Variance Blocking Variable (Subject): Subject

Group: Treatment

Now let’s see.


                                             PE       90% CI     CVtot. CV(R) CV(T)

orig. data set R 3.0.1     var.equal=TRUE   83.66  63.51 110.19

                           var-equal=FALSE  83.66  63.49 110.22

               Phoenix 6.3 standard coding  83.66  63.51 110.19  40.86

                           modified coding  83.66  63.49 110.22  40.86  37.14 44.35

mod. data set  R 3.0.1     var.equal=TRUE  124.35  81.21 190.41

                           var-equal=FALSE 124.35  76.36 202.51

               Phoenix 6.3 standard coding 124.35  81.21 190.41  60.54

                           modified coding 124.35  76.36 202.51  60.54  37.14 87.04

As a side effect Phoenix spits out the treatments’ total variances in Var(Period*Treatment*Subject)_11 (R) and Var(Period*Treatment*Subject)_12 (T). In the balanced case the mean of these variances equals the total (pooled) total variance of the standard set-up. ~~Have to figure out yet how to~~ weight them properly in the imbalanced case.*

In the next release of Phoenix (~1^st half 2014), the dummy variable will not be needed any more (i.e., can be left empty).

<div lang="de" xml:lang="de">

Was lange währt, wird endlich gut.

</div>

s²_tot.=(s²_Rν_R+s²_Tν_T)/(ν_R+ν_T), where ν_(R,T) = n_(R,T) – 1.
Quick & dirty R:
round(100*sqrt(exp(var(logAUC[treatment=="R"]))-1), 2) round(100*sqrt(exp(var(logAUC[treatment=="T"]))-1), 2) round(100*sqrt(exp((var(logAUC[treatment=="R"])*(length(logAUC[treatment=="R"])-1)+ var(logAUC[treatment=="T"])*(length(logAUC[treatment=="T"])-1))/ (length(logAUC[treatment=="R"])+length(logAUC[treatment=="T"])-2))-1), 2)

—
Dif-tor heh smusma 🖖🏼 Довге життя Україна!
Helmut Schütz

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

ElMaestro
★★★

Denmark,
2013-07-04 23:16
(4309 d 17:44 ago)

@ Helmut
Posting: # 10939
Views: 11,917

Parallel design [unequal variances] in Phoenix/WinNonlin

Post reply

Hi Helmut and all,

this looks to me like an ugly work-around rather than a solution.
I am not a WNL user but I'd like to understand why that type of coding fixes it. I believe WNL will fit this model via REML, and on basis of the syntax I imagine this means a covariance matrix with two squared sigmas on the diagonal and zero elsewhere. But to get through to the 'correct' CI there must be a step for derivation of the denominator DF which happens to coincide with the DF you get via the Welch correction when just doing ordinary unpaired t-test with unequal variances. Perhaps that's similar to the SATT option in SAS (which I also don't use). Could be checked?

So, what is really going on within WNL when that weird coding principle is applied - Why does it give the seemingly right result? Anyone knows?

—
Pass or fail!
ElMaestro

Helmut
★★★

Vienna, Austria,
2013-07-05 00:13
(4309 d 16:47 ago)

@ ElMaestro
Posting: # 10940
Views: 11,868

Workaround

Post reply

Hi ElMaestro,

❝ this looks to me like an ugly work-around rather than a solution.

Right. Solution was my term. ;-)

Linda wrote:

Set up the Repeated subtab of the Variance Structure tab as:
Repeated Specification: <blank> (see below)
Variance Blocking Variables (Subject): Subject
Group: Formulation
Type: Variance Components
However, in the current version of Phoenix, if the Repeated Specification is blank, Phoenix 1.3 will ignore the rest of this tab, so a workaround must be used. […] After completion of these steps, execution of Bioequivalence will produce the same values as SAS.

❝ […] I'd like to understand why that type of coding fixes it.

So do I.

❝ I believe WNL will fit this model via REML, …

REML is the only language WinNonlin speaks.

❝ … and on basis of the syntax I imagine this means a covariance matrix with two squared sigmas on the diagonal and zero elsewhere.

You are the Matrazenkaiser. Reasonable assumption.

❝ But to get through to the 'correct' CI there must be a step for derivation of the denominator DF which happens to coincide with the DF you get via the Welch correction when just doing ordinary unpaired t-test with unequal variances. Perhaps that's similar to the SATT option in SAS (which I also don't use). Could be checked?

Sure. Satterthwaite is WinNonlin’s standard setup. I’m not in the office so can’t check what will happen with the other option ‘Residual’.

❝ So, what is really going on within WNL when that weird coding principle is applied - Why does it give the seemingly right result? Anyone knows?

Maybe. Too many cooks spoil the broth. Pharsight’s statisticians ≠ coders.

—
Dif-tor heh smusma 🖖🏼 Довге життя Україна!
Helmut Schütz

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

ElMaestro
★★★

Denmark,
2013-07-05 00:37
(4309 d 16:23 ago)

@ Helmut
Posting: # 10941
Views: 11,876

3-trt study?

Post reply

Hi Helmut,

thanks for the comments.
In perspective, could this methodology then be extended to a 3-treatment parallel study (e.g. US Ref, EU Ref, Test)? I am thinking we'd just need another sigma squared on the diagonal and then be good to go? If that's the case then the WNL approach –however counterintuitive and retarded it might seem– could provide a solution to a problem which cannot be resolved by the standard t-test approach?
Oh and by the way, I think this would imply that the denominator DFs would have to be calculated from treatment A, B as well as C, when a 90% CI for A/B is constructed - that's perhaps also worthy of a heated debate further down the road...?

:pirate:

—
Pass or fail!
ElMaestro

Helmut
★★★

Vienna, Austria,
2013-07-05 01:08
(4309 d 15:52 ago)

@ ElMaestro
Posting: # 10942
Views: 11,958

؟ῳ ♠♣♥♦♪♫ѳ٭

Post reply

Hi ElMaestro,

❝ In perspective, could this methodology then be extended to a 3-treatment parallel study (e.g. US Ref, EU Ref, Test)?

Intuitively – if we walk on EMA’s stony ground of ‘leave one out’ – I would say so.* Duno what will happen if we keep all treatments in the model.

❝ I am thinking we'd just need another sigma squared on the diagonal and then be good to go?

I feel tempted to give you all statisticians’ standard answer on any arbitrary question: “Positively maybe!”

❝ If that's the case then the WNL approach –however counterintuitive and retarded it might seem– could provide a solution to a problem which cannot be resolved by the standard t-test approach?

See above. Any ideas how to code that in R without setting the Silly-o-Meter on fire?

❝ Oh and by the way, I think this would imply that the denominator DFs would have to be calculated from treatment A, B as well as C, when a 90% CI for A/B is constructed - that's perhaps also worthy of a heated debate further down the road...?

Sure. We are moving in a circle. Remember the sim’s Detlew performed in Dec 2011?

If I erred, I am prepared to accept discipline up to level three.
Punishment of a higher order will be considered unjustified
and rejected with prejudice. David Brin (Startide Rising, p30, 1983)

—
Dif-tor heh smusma 🖖🏼 Довге життя Україна!
Helmut Schütz

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

Helmut
★★★

Vienna, Austria,
2013-07-05 17:13
(4308 d 23:47 ago)

@ ElMaestro
Posting: # 10945
Views: 12,466

Satterthwaite’s df in Phoenix/WinNonlin

Post reply

Hi ElMaestro,

❝ […] there must be a step for derivation of the denominator DF which happens to coincide with the DF you get via the Welch correction when just doing ordinary unpaired t-test with unequal variances.

Given the two variances and sample sizes PHX/WNL’s degrees of freedom agree with R (21.4306 for the full data set and 10.7555 for the modified one); checked to 14 significant digits.

According to Phoenix’ User’s Guide

The linear model used in linear mixed effects modeling is:

Y~N(Xb,V)

where X is a matrix of fixed known constants, b is an unknown vector of constants, and V is the variance matrix dependent on other unknown parameters. Wald (1943)¹ developed a method for testing hypotheses in a rather general class of models. To test the hypothesis H₀: Lb = 0, the Wald statistic is:

(Lβ)^T[L([X^TV^–1X)^–1L^T]^–1(Lβ)

where β and V are estimators of β and V. L must be a matrix such that each row is in the row space of X. The Wald statistic is asymptotically distributed under the null hypothesis as a χ² random variable with q = rank(L) degrees of freedom. In common variance component models with balanced data, the Wald statistic has an exact F-distribution. Using the χ² approximation results in a test procedure that is too liberal. LinMix uses a method based on analysis techniques described by Allen and Cady (1982)² for finding an F-distribution to approximate the distribution of the Wald statistic.
Giesbrecht and Burns (1985)³ suggested a procedure to determine denominator degrees of freedom when there is one numerator degree of freedom. Their methodology applies to variance component models with unbalanced data. Allen derived an extension of their technique to the general mixed model.

… followed by two pages in matrix-notation I can’t comprehend (aka hazelnut brain).

Let’s have a look at the modified data set (full precision).
s²₁ 0.129240461704103, n₁ 12
s²₂ 0.563916088422299, n₂ 9
If I trust Wikipedia
(s²₁+s²₂)²/[s⁴₁/(n₁–1)+s⁴₂/(n₂–1)] = 11.64240178 (≠ 10.75546079) :confused:

Aha!
(s²₁/n₁+s²₂/n₂)²/{s⁴₁/[n₁²(n₁–1)]+s⁴₂/[n₂²(n₂–1)]} = 10.75546079

[image]

Funny enough the radio buttons above turn out to be a gimmick without any effect. Whatever you choose – the dfs are identical.

Strange:

[image]

Wald A. Tests of statistical hypotheses concerning several parameters when the number of observations is large. Transactions of the American Mathematical Society 1943;54(3):426–82.*
Allen DM, Cady FB.[/b] Analyzing Experimental Data By Regression. Belmont: Lifetime Learning Publications; 1982.
Giesbrecht FG, Burns JC. Two-Stage Analysis Based on a Mixed Model: Large-Sample Asymptotic Theory and Small-Sample Simulation Results. Biometrics. 1985;41(2):477–86.

Wow, 57 pages! This must be serious and important stuff since in comparison Mr Einstein needed only 31 pages to describe special relativity (in German, which is more verbose than English)…
Einstein A. Zur Elektrodynamik bewegter Körper. Annalen der Physik. 1905;4(17–5): 891–921.

—
Dif-tor heh smusma 🖖🏼 Довге життя Україна!
Helmut Schütz

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

ElMaestro
★★★

Denmark,
2013-07-05 19:08
(4308 d 21:52 ago)

@ Helmut
Posting: # 10946
Views: 11,920

Satterthwaite’s df in Phoenix/WinNonlin

Post reply

Hi Helmut,

❝ If I trust Wikipedia

❝ (s²₁+s²₂)²/[s⁴₁/(n₁–1)+s⁴₂/(n₂–1)] = 11.64240178 (≠ 10.75546079) :confused:

10.76 appears correct.
Quickie:

WelchDF=function(V1, V2, N1, N2)

{

 Num = ((V1/N1)+(V2/N2))^2

 Denom = (V1/N1)^2/(N1-1) + (V2/N2)^2/(N2-1) 

 return(Num/Denom)

}

But I have to admit I am still one big question mark as to why all this works correctly. All the matrix explanation and stuff is rather gibberish to me. I looked up what SAS writes about DDFM=SATT and I felt like hitting myself with a hammer after reading two lines of it.

—
Pass or fail!
ElMaestro

Helmut
★★★

Vienna, Austria,
2013-07-05 19:15
(4308 d 21:46 ago)

@ ElMaestro
Posting: # 10947
Views: 11,885

Wikipedia, the free encyclopedia that anyone can edit.

Post reply

Hi ElMaestro,

❝ ❝ If I trust Wikipedia

❝ ❝ (s²₁+s²₂)²/[s⁴₁/(n₁–1)+s⁴₂/(n₂–1)] = 11.64240178 (≠ 10.75546079) :confused:

WP is a strange piece of cake. The formula was correct on 2013-04-16,

[image]

an anonymous editor from Austin, Texas screwed up on 2013-06-27,

[image]

and corrected today by my humble self.

[image]

I love audit-trails.

❝ 10.76 appears correct.

In the meantime I have edited my post…

❝ Quickie:

Nice.

❝ But I have to admit I am still one big question mark as to why all this works correctly. All the matrix explanation and stuff is rather gibberish to me. I looked up what SAS writes about DDFM=SATT and I felt like hitting myself with a hammer after reading two lines of it.

You are not alone.

—
Dif-tor heh smusma 🖖🏼 Довге життя Україна!
Helmut Schütz

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

ElMaestro
★★★

Denmark,
2013-07-05 20:25
(4308 d 20:35 ago)

@ Helmut
Posting: # 10948
Views: 11,907

Wikipedia, the free encyclopedia that anyone can edit.

Post reply

Hi Helmut,

at this point allow me to sum it up so far:

Noone knows why, but if you are an egg-head with an IQ of 6 billion and qualified for a professorship at MIT and you chose to play around with 9 incomprehensible matrices from an obscure mixed model whose specification in WinNonlin is completely counterintuitive then you may be able to achieve exactly the same results as any lay person like you and I can derive from using 4 seconds on googling "Welch-Satterthwaite".

Science is really making progress these days.

Pharsight are you reading this? (Come on Simon, we know you are!)

—
Pass or fail!
ElMaestro

Helmut
★★★

Vienna, Austria,
2013-07-05 20:40
(4308 d 20:20 ago)

@ ElMaestro
Posting: # 10949
Views: 12,097

Improbable

Post reply

Hi ElMaestro,

❝ […] if you are an egg-head with an IQ of 6 billion and qualified for a professorship at MIT…

At least chances for the former are somewhat low:

pnorm(6e9, mean=100, sd=15, lower.tail=FALSE)

❝ Science is really making progress these days.

Yep.

❝ Pharsight are you reading this? (Come on Simon, we know you are!)

BTW, to quote Senn and Grieve*

(2) Why would an investigator plan a trial with the object of proving equality of two formulations if the variances were believed different?
(4) Welch’s test is not universally accepted as being a valid approach to the problem of statistical inference for differences in means in the presence of heteroscedasticity.

Oh boy!

Senn S, Grieve AP. A Comment on Optimal Allocations for Bioequivalence Studies. Biometrics. 1999;55(4):1314–5.doi 10.1111/j.0006-341X.1999.01314.x

—
Dif-tor heh smusma 🖖🏼 Довге життя Україна!
Helmut Schütz

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