Splus codes applied to a simulated data set for tests
  in "Assessing Antiviral Potency of Anti-HIV Therapies in Vivo
         by Comparing Viral Decay Rates in Viral Dynamic Models",
                    by A. Adam Ding and Hulin Wu

Notes:

1. The data were randomly generated similar to the schedules
   of the real data set analyzed in Ding and Wu (2000, Biostatistics).

2. A total of 32 patients are contained in this study, indicated by
   the ID variable. The Day variable reflects the time for viral
   load measurement time after initiation of antiviral treatment.
   The RNA variable contains the viral load measurement in unit of
   copies per ml blood. The patients are separated into two
   treatment groups (Group 1 and Group 2) with 16 patients each,
   and is identified in the Group variable.

3. Detection limit of the viral load (HIV RNA copies) assay
   is 400 copies per ml blood. If it is below detectable, the data
   is reflected as 400. We imputed 200 for these data points in the
   following data analysis. If more than one measurement below
   detectable level for an individual, we just impute the first
   measurement and exclude the rest of them in our analysis (otherwise
   it may result in misleading dynamic patterns).
 

4. Due to convergence problems, multiple starting values
   (or global search algorithms) should be used to fit the
   models. The starting point in the following Splus code was
   chosen after comparing the results using many different
   starting values.
 

######################### Beginning of Splus Code ########################
# Input Data: assume that the data at the end is stored in an ASCII file
# named "data".
workd<-read.table(file="data",header=T,row.names=NULL)

# Impute below detectable viral measurement by half of the detection limit.
workd$RNA[workd$RNA==400]<-200

# For patient number 27, the viral measurements are below detectable level on
# both 14th and 28th day. So the latter point is excluded from the analysis.
workd<-workd[(workd$ID!=27 | workd$Day!=28),]

# Define a function representing the 2 phrase decay model (1.2)
f.twoe<-function(p1, d1, p2, d2, X)
{
        P1 <- exp(p1)
        P2 <- exp(p2)
        P1 * exp( - d1 * X) + P2 * exp( - d2 * X)
}

# Fit the PMEM
# We used a fixed starting values of p1=12, d1=0.65, p2=8 and d2=0.03
# on this data set.
# For other data sets, several starting values should be tried and
# then the best fit chosen.
tmp.fit<-nlme(Y ~ log10(f.twoe(p1, d1, p2, d2, X)),
              fixed  = list(p1~.,  d1~., p2~.,  d2~.),
              random = list(p1~., d1~., p2~., d2~.),
              cluster = ~ Z,
              data = data.frame(Y=log10(workd$RNA),
                     X=workd$Day, Z=workd$ID),
              start = list(fixed= c(12, 0.65,  8, 0.03)))

# Get the EBEs for the decay rates.
tmp.rate<-t(t(tmp.fit$coe$random)+tmp.fit$coe$fixed)[,c("d1","d2")]

# Separate the EBEs into the groups to prepare for hypothesis testing
patid<-sort(unique(workd$ID))
np<-length(patid)
group<-(1:np)
for (i in 1:np) {group[i]<-workd$Group[workd$ID==patid[i]][1]}

 np1<-sum(group==1)
 np2<-sum(group==2)
 d1.grp1<-tmp.rate[group==1,"d1"]
 d1.grp2<-tmp.rate[group==2,"d1"]
 d2.grp1<-tmp.rate[group==1,"d2"]
 d2.grp2<-tmp.rate[group==2,"d2"]

# Conduct EBE-based tests:
# t-test and rank-sum test for AH1
# t-test and rank-sum test for AH2
# T^2-test (MANOVA) and O'Brien rank-sum test for AH3

# AH1. Testing one-sided alternative here: d1 is bigger in Group 2.
t.test(d1.grp1,d1.grp2,alternative="less")$p.value
wilcox.test(d1.grp1,d1.grp2,alternative="less")$p.value
# Resulting p-values: t.test 0.049997, wilcox.test 0.09175338

# AH2. Testing one-sided alternative here: d2 is bigger in Group 2.
t.test(d2.grp1,d2.grp2,alternative="less")$p.value
wilcox.test(d2.grp1,d2.grp2,alternative="less")$p.value
# Resulting p-values: t.test 0.4378882, wilcox.test 0.4482666

# AH3. Testing two-sided alternative: d1 and/or d2 are different.
# MANOVA
summary(manova(rbind(cbind(d1.grp1,d2.grp1),cbind(d1.grp2,d2.grp2)) ~ group))$Stats[,1,5]
# Resulting p-value 0.2614775

# O'Brien test
 tmp<-rank(c(d1.grp1,d1.grp2))
 rd1.grp1<-tmp[(1:np1)]
 rd1.grp2<-tmp[-(1:np1)]
 tmp<-rank(c(d2.grp1,d2.grp2))
 rd2.grp1<-tmp[(1:np1)]
 rd2.grp2<-tmp[-(1:np1)]
 rd.grp1<-rd1.grp1+rd2.grp1
 rd.grp2<-rd1.grp2+rd2.grp2
wilcox.test(rd.grp1,rd.grp2)$p.value
# Resulting p-value 0.2344624

###################### End of Splus Code #####################################
 
 
 

ID   Day    RNA   Group
 1  0.000    67324 1
 1  0.125   171810 1
 1  0.333    39576 1
 1  1.000    56412 1
 1  3.000     5688 1
 1  7.000     1956 1
 1 14.000     1379 1
 1 28.000      400 1
 2  0.000   569633 1
 2  0.125   248890 1
 2  0.333   425369 1
 2  1.000   119297 1
 2  3.000    94723 1
 2  7.000     9045 1
 2 14.000      867 1
 2 28.000      805 1
 3  0.000   398899 1
 3  0.125   219505 1
 3  0.333   178366 1
 3  1.000    73744 1
 3  3.000    35458 1
 3  7.000    51040 1
 3 14.000     1701 1
 3 28.000     1078 1
 4  0.000   610484 1
 4  0.125   323648 1
 4  0.333   271977 1
 4  1.000   598412 1
 4  3.000    72110 1
 4  7.000     7902 1
 4 14.000     4696 1
 4 28.000      761 1
 5  0.000   253495 1
 5  0.125   432427 1
 5  0.333   367006 1
 5  1.000    50973 1
 5  3.000    41538 1
 5  7.000     9444 1
 5 14.000      991 1
 5 28.000      711 1
 6  0.000  1541533 1
 6  0.125   860756 1
 6  0.333   593412 1
 6  1.000   339032 1
 6  3.000   181938 1
 6  7.000    76764 1
 6 14.000    19341 1
 6 28.000     5218 1
 7  0.000    23678 1
 7  0.125    70048 1
 7  0.333   138385 1
 7  1.000    32492 1
 7  3.000    11417 1
 7  7.000     1931 1
 7 14.000      421 1
 7 28.000      400 1
 8  0.000   744825 1
 8  0.125   521558 1
 8  0.333   325766 1
 8  1.000   356948 1
 8  3.000    52540 1
 8  7.000    11829 1
 8 14.000     3254 1
 8 28.000    12596 1
 9  0.000  1252328 1
 9  0.125   261470 1
 9  0.333   320139 1
 9  1.000   142677 1
 9  3.000    14074 1
 9  7.000    12346 1
 9 14.000     8613 1
 9 28.000     3363 1
10  0.000   329416 1
10  0.125   168726 1
10  0.333   291034 1
10  1.000   111688 1
10  3.000    70522 1
10  7.000    14427 1
10 14.000     1548 1
10 28.000     1282 1
11  0.000   157635 1
11  0.125    85405 1
11  0.333    97043 1
11  1.000    66889 1
11  3.000     6094 1
11  7.000     2765 1
11 14.000      995 1
11 28.000      400 1
12  0.000   463754 1
12  0.125  1111460 1
12  0.333   461062 1
12  1.000   340430 1
12  3.000   166770 1
12  7.000    57106 1
12 14.000     9993 1
12 28.000     1856 1
13  0.000   225887 1
13  0.125   177991 1
13  0.333   109232 1
13  1.000   249347 1
13  3.000    42956 1
13  7.000     2712 1
13 14.000     1344 1
13 28.000      400 1
14  0.000  5354718 1
14  0.125 10156869 1
14  0.333  8465124 1
14  1.000  6276483 1
14  3.000  1759661 1
14  7.000   372656 1
14 14.000    67606 1
14 28.000    33587 1
15  0.000   162043 1
15  0.125    47688 1
15  0.333   114639 1
15  1.000    67453 1
15  3.000    16189 1
15  7.000     4507 1
15 14.000     1269 1
15 28.000      441 1
16  0.000   466821 1
16  0.125   323978 1
16  0.333   288308 1
16  1.000   134974 1
16  3.000    33357 1
16  7.000     6074 1
16 14.000     2573 1
16 28.000     3130 1
17  0.000   766694 2
17  0.125   595493 2
17  0.333   287100 2
17  1.000   435330 2
17  3.000    28789 2
17  7.000     2615 2
17 14.000     2145 2
17 28.000     1867 2
18  0.000   154188 2
18  0.125   109839 2
18  0.333   200805 2
18  1.000   151538 2
18  3.000    37298 2
18  7.000     2224 2
18 14.000      490 2
18 28.000      616 2
19  0.000    70841 2
19  0.125   252774 2
19  0.333   192880 2
19  1.000   151659 2
19  3.000    43354 2
19  7.000     6732 2
19 14.000     1738 2
19 28.000     1377 2
20  0.000   158011 2
20  0.125   183884 2
20  0.333   114733 2
20  1.000   139015 2
20  3.000    29253 2
20  7.000     4070 2
20 14.000     1868 2
20 28.000      672 2
21  0.000   934817 2
21  0.125  1352673 2
21  0.333  1192666 2
21  1.000   142216 2
21  3.000    55779 2
21  7.000    32402 2
21 14.000     3525 2
21 28.000     1039 2
22  0.000   227665 2
22  0.125   123610 2
22  0.333    98786 2
22  1.000    60500 2
22  3.000    13351 2
22  7.000     1044 2
22 14.000      691 2
22 28.000      400 2
23  0.000  2961627 2
23  0.125   554609 2
23  0.333   902410 2
23  1.000   676630 2
23  3.000   164742 2
23  7.000    45623 2
23 14.000    13675 2
23 28.000     8694 2
24  0.000   421949 2
24  0.125   136907 2
24  0.333   732840 2
24  1.000   133798 2
24  3.000    16982 2
24  7.000     9185 2
24 14.000     3639 2
24 28.000     1401 2
25  0.000   926362 2
25  0.125  1488839 2
25  0.333   364706 2
25  1.000   663691 2
25  3.000   100743 2
25  7.000    22729 2
25 14.000     7972 2
25 28.000     1218 2
26  0.000   302650 2
26  0.125   359009 2
26  0.333   599914 2
26  1.000   323553 2
26  3.000   116214 2
26  7.000     2533 2
26 14.000     6498 2
26 28.000     2381 2
27  0.000    59804 2
27  0.125    22617 2
27  0.333    22926 2
27  1.000    22569 2
27  3.000     3523 2
27  7.000     1507 2
27 14.000      400 2
27 28.000      400 2
28  0.000   284082 2
28  0.125   101536 2
28  0.333   134319 2
28  1.000    57690 2
28  3.000    20888 2
28  7.000     1340 2
28 14.000      634 2
28 28.000      524 2
29  0.000    50026 2
29  0.125   146659 2
29  0.333    64786 2
29  1.000    43144 2
29  3.000    12406 2
29  7.000     2341 2
29 14.000      640 2
29 28.000      486 2
30  0.000    83640 2
30  0.125   111536 2
30  0.333    84267 2
30  1.000    91913 2
30  3.000    24772 2
30  7.000     6968 2
30 14.000     1028 2
30 28.000      911 2
31  0.000   345629 2
31  0.125   233055 2
31  0.333   287857 2
31  1.000    52375 2
31  3.000    38088 2
31  7.000     4804 2
31 14.000     3548 2
31 28.000     2749 2
32  0.000   867389 2
32  0.125   523124 2
32  0.333   244840 2
32  1.000   174122 2
32  3.000    44033 2
32  7.000     7083 2
32 14.000     1445 2
32 28.000      421 2