Bioequivalence and Bioavailability Forum

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 (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)