Vignette Info

A detailed description of the botulinum study is provided on p 369-370 in the text. This is a multivariate case-control study with missing data and repeated measures. The aim of the study is to assess the effect of botox, an injection used to treat myofascial pain. Several variables are measured as outcomes; they are described in the text. The null hypothesis is that the botox treatment had no effect, and the alternative hypothesis is that patients treated with botox had lower/higher scores than those in the placebo group (one-sided alternatives). To account for baseline measurements, we consider the difference between observations at time $t=1,2,3$ and the observations at time $t=0$ rather than the raw measurements.

To begin, we read the data, collect informations about samples (individual), time (time), an indicator for whether the patient received botox (botox) or not (sol_fisio) (treatment), and separate them from the outcomes (covariates).

library(permuter)
library(knitr)
data(botulinum)
set.seed(100)

individual <- botulinum[1:20 , "paz"]
treatment <- botulinum[1:20 , "farmaco"]
time <- botulinum[ , "Time"]
covariates <- botulinum[, -(1:4)]
n <- length(unique(individual))
p <- ncol(covariates)

We abstract the testing problem: a $V$-dimensional non-degenerate variable is observed at $k$ different times on $n_0$ and $n_1$ units in two experimental situations, corresponding to two levels of a symbolic treatment. In our study, $V=24, k=4, n_0 = n_1 = 10$, and $n = n_0+n_1 = 20$. In this longitudinal study there is a main obstacle to applying classic parametric tests: the number of observed variables at different time points is much higher than the number of subjects ($V \cdot k \gg n$). Furthermore, since all variables may be informative for differentiating two groups, the NPC approach properly applies when analyzing these data. Classic parametric tests or even rank tests in such situations may fail to take into account the dependence structure across variables and time points.

The whole data set is denoted by: $$ \begin{aligned} \mathbf{X} &= { X_{hji}(t), t=1, \dots, k, i=1,\dots,n_j, j=0,1, h=1,\dots,V} \ &= {\mathbf{X}{hji}, i=1,\dots,n_j, j=0,1, h=1,\dots,V} \end{aligned} $$ where $\mathbf{X}{hji} = { X_{hji}(t), t = 1,\dots,k}$.

In order to account for different baseline observations, we consider the $k - 1$ $V$-dimensional differences $D_{hji}(t) = X_{hji}(t) - X_{hji}(0), t = 1, 2, 3$, $i = 1,\dots,n_j$, $j = 0,1$, $h = 1,\dots,V$. Hence the hypothesis testing problem related to the $h$th variable may be formalized as

$$H_{0h}: \left\lbrace \bigcap_{t=2}^k \left[ D_{h0}(t) \,{\buildrel d \over =}\, D_{h1}(t)\right] \right\rbrace = \left\lbrace \cap_{t=2}^k H_{0ht}\right\rbrace, h=1,\dots,V$$ against the alternative: $$H_{1h}: \left\lbrace \bigcup_{t=2}^k H_{1ht}\right\rbrace, h=1,\dots,V,$$

where $H_{1ht}: \left[ D_{h0}(t) \,{\buildrel d \over >}\,D_{h1}(t)\right]$ or $\left[ D_{h0}(t) \,{\buildrel d \over <}\, D_{h1}(t)\right]$ according to which kind of stochastic dominance is of interest for the $h$th variable. The alternative hypothesis is that patients treated with the botulinum toxin had lower values than those treated with the placebo ("less than" alternative) for the first nineteen variables, and that patients treated with botulinum toxin had higher values than those in the placebo group for the variables Mas, Maf, Mp, Mld, and Mls ("greater than" alternative). We encode these alternative hypotheses in the vector alternatives.

colnames(covariates)
greater_than <- c("Mas", "Maf", "Mp", "Mld", "Mls")
alternatives <- ifelse(colnames(covariates) %in% greater_than, "greater", "less")

First Analysis

We create a list called D, where each element of the list is an \Sexpr{n}$\times$\Sexpr{p} matrix of measurements corresponding to a different time point. We compute the differences $D_{ijk}$ with respect to time zero observations and remove the vector of observations at $t_0$ (D[[1]]). The original data is ordered so that the first \Sexpr{n} rows correspond to time 0, the next \Sexpr{n} rows correspond to time 1, and so forth. Furthermore, the individuals are stored in the same order within each of these \Sexpr{n}-row sections. This facilitates the data processing.

D <- list()
for (t in 0:3) {
    D[[t+1]] <- covariates[time == t, ]
}

# Take differences, Tj - T0
for (t in 2:4) {
    D[[t]] <- D[[t]] - D[[1]]
}
D <- D[-1]

We first assess the effect of botulinum by performing a two-sample test with repeated measures. We use a modified test statistic, essentially a weighted sum of the difference in means between treatment groups over time. The weights depend on the number of non-missing observations. Since the distribution of non-missing observations varies between variables, we must consider one variable at time. The vector o is the indicator function (missing/not missing) applied to each variable and the vector nu contains the number of non-missing observations in the two samples. The test statistic for the $j$th variable is: $T_j = \sum^3_{k=1}[\omega_j D_{1jk} - \omega_j D_{0jk} ]$, where $\omega_j = \sqrt{\nu_{0j}/\nu_{1j}}$.

count_nas <- function(covariate, treatment){
  tr <- unique(treatment)
  nu <- rep(0, 2)
  observed <- !is.na(covariate)
  nu[1] <- sum(observed & treatment == tr[1])
  nu[2] <- sum(observed & treatment == tr[2])
  return(nu)
}


compute_test_statistic <- function(covariate_list, treatment){
  # covariate_list: a list where each element is a matrix/dataframe 
  # containing the measurements at a particular time point.
  # The dimensions of the matrices must be the same and observations
  # should be ordered the same in each matrix.
  # treatment: a vector indicating which of two treatments each 
  # individual received. Same order as the individuals in each matrix.
  #
  # Returns the test statistic T_j for each variable (column).

  timepoints <- length(covariate_list)
  p <- ncol(covariate_list[[1]])

  treat_levels <- unique(treatment)
  tst <- rep(0, p)

  for (t in seq_len(timepoints)) {
    group1 <- covariate_list[[t]][treatment == treat_levels[1], ]
    group0 <- covariate_list[[t]][treatment == treat_levels[2], ]
    for (j in seq_len(p)) {
      nu <- count_nas(covariate_list[[t]][, j], treatment)
      omega <- sqrt(nu[2]/nu[1])
      tst[j] <- tst[j] + omega * (mean(group1[, j], na.rm = TRUE) - mean(group0[, j], na.rm = TRUE))
    }
  }
  return(tst)
}

Under the null hypothesis of no effect of botox, treatment is essentially a meaningless label and so the treatment assignments are exchangeable. In order to obtain the null distribution, we permute observations by randomly exchanging the rows of the observations at each time point, applying the same permutation to each time point. This ensures that the inner dependencies among variables and due to the repeated mesures are maintained. Note, this is equivalent to simply permuting the treatment assignments, but conceptually emphasizes the need to preserve dependencies.

B <- 1000

tst <- compute_test_statistic(D, treatment)

distr <- t(replicate(B, {
  D_star <- permute_rows(D)
  compute_test_statistic(D_star, treatment)
}))

We obtain the raw p-values using t2p and then adjust for multiplicity using fwe_minp with Tippett's combining function.

partial_pvalues <- sapply(1:p, function(j) t2p(tst[j], distr[, j], alternatives[j]))
adjusted_pvalues <- fwe_minp(partial_pvalues, distr, combine = "tippett")

res <- data.frame(colnames(covariates), partial_pvalues, adjusted_pvalues)
colnames(res) <- c("Variable", "Raw p-values", "Adjusted p-values")
kable(res)

Second Analysis

In contrast to the first analysis, we focus on only on the individuals who received botox and look at the differences between measurements at each time point.

NB: there is no description of this analysis anywhere. Not sure what is being tested.

Time <- botulinum$Time+1
paz <- (treatment=="botox")
D <- array(0, dim = c(n, p, 4))

for (t in 1:4) {
    D[, , t] <- as.matrix(covariates[Time == t, ])
}

### differences T_j-T_{j-1}

for (t in 2:4) {
    D[, , t] <- D[, , t] - D[, , (t - 1)]
}

D <- D[, , -1]



DB <- D[paz == 1, , ]
n <- dim(DB)[1]

############################################################################## NEW Take differences over time for individuals who received botox
D_new <- list()
for (t in 0:3) {
    D_new[[t+1]] <- covariates[time == t & treatment == "botox", ]
}

# Take differences, T_j-T_{j-1}
for (t in 2:4) {
    D_new[[t]] <- D_new[[t]] - D_new[[t - 1]]
}
D_botox <- D_new[-1]
n <- nrow(D_botox[[1]])
############################################################################# END NEW

The test statistic we use is $T_j = \sum_{k=1}^3 \left( S_{jk} \sqrt{ \frac{ \sum_{i\neq k} \nu_i}{\nu_k}} - (\sum_{i=1}^3 S_{ji}-S_{jk})\sqrt{ \frac{ \sum_{i\neq k} \nu_i}{\nu_k}}\right)^2$ where $S_{jk}$ is the sum of non-missing differences for the $j$th variable at the $k$th difference in times and $\nu_k$ is the number of missing values for the $j$th variable at the $k$th difference in times.

# Old method
B <- 1000
T <- array(0, dim = c((B + 1), p))
n <- 10

for (j in 1:p) {

    Y <- DB[, j, ]
    O <- apply(Y, 2, function(x) {
        ifelse(is.na(x), 0, 1)
    })
    Y <- ifelse(is.na(DB[, j, ]), 0, DB[, j, ])

    nu <- apply(O, 2, sum)
    S <- apply(Y, 2, sum)

    for (t in 1:3) {

        T[1, j] <- T[1, j] + (S[t] * sqrt(sum(nu[-t])/nu[t]) - (sum(S) - S[t]) * 
            sqrt(sum(nu[-t])/nu[t]))^2

    }
}  #end p



for (bb in 2:(B + 1)) {

    D.star <- apply(DB, c(1, 2), sample)
    D.star <- aperm(D.star, c(2, 3, 1))
    for (j in 1:p) {

        Y <- D.star[, j, ]
        O <- apply(Y, 2, function(x) {
            ifelse(is.na(x), 0, 1)
        })
        Y <- ifelse(is.na(D.star[, j, ]), 0, D.star[, j, ])

        nu <- apply(O, 2, sum)
        S <- apply(Y, 2, sum)
        for (t in 1:3) {

            T[bb, j] <- T[bb, j] + (S[t] * sqrt(sum(nu[-t])/nu[t]) - (sum(S) - S[t]) * 
                sqrt(sum(nu[-t])/nu[t]))^2

        }
    }  #end p

}


P <- t2p_old(T)

P[1, ]

Not sure if this is right -- the p-values don't match at the end, but I can't find an issue with the math

############################################################################# NEW
tst <- rep(0, p)
for (j in 1:p) {
    Y <- sapply(D_botox, function(x) x[, j])
    nu <- apply(Y, 2, function(x) sum(!is.na(x)))
    S <- apply(Y, 2, sum, na.rm = TRUE)

    for (t in 1:3) {
        tst[j] <- tst[j] + (S[t] * sqrt(sum(nu[-t])/nu[t]) - (sum(S) - S[t]) * sqrt(sum(nu[-t])/nu[t]))^2
        # tst[j] <- tst[j] + (sum(nu[-t])/nu[t])*(2*S[t]-sum(S))^2 # Equivalent, fewer
        # calculations
    }
}

time_ints <- rep(1:3, each = n)
D_temp <- do.call(rbind, D_botox)  # combine the 3 matrices into one
distr <- matrix(0, nrow = B, ncol = p)
for (bb in 1:B) {
    # time_ints <- sample(time_ints) # permute the times: under the null, the
    # distributions for each are equal
    D_star <- apply(D_temp, 2, sample)  # permute within columns, break the correlation within columns
    D_star <- split(data.frame(D_star), time_ints)  # break up the combined data into 3 random matrices of equal size, preserving rows
    for (j in 1:p) {
        Y <- sapply(D_star, function(x) x[, j])
        nu <- apply(Y, 2, function(x) sum(!is.na(x)))
        S <- apply(Y, 2, sum, na.rm = TRUE)
        for (t in 1:3) {
            distr[bb, j] <- distr[bb, j] + (S[t] * sqrt(sum(nu[-t])/nu[t]) - (sum(S) - 
                S[t]) * sqrt(sum(nu[-t])/nu[t]))^2  #+ (sum(nu[-t])/nu[t])*(2*S[t]-sum(S))^2 # Equivalent, fewer calculations       
        }
    }
}

raw_pvalues <- sapply(1:p, function(j) t2p(tst[j], distr[, j], alternative = "greater"))

############################################################################# END NEW

kable(cbind(raw_pvalues, `old method p-values` = P[1, ]), row.names = FALSE)


statlab/permuter documentation built on May 30, 2019, 9:45 a.m.