pairedCompMean: Permutation Test for the Mean with Paired Data

In bozenne/butils: Various useful functions

Permutation Test for the Mean with Paired Data

Description

Permutation test for comparing mean between within cluster with efficient adjustment for multiple comparisons

Usage

pairedCompMean(
  data,
  col.treatment,
  col.strata,
  col.value,
  col.cluster,
  n.perm = 10000,
  method.adj = "single-step",
  seed = NULL,
  trace = TRUE
)

Arguments

data	[data.frame] dataset with in row the samples.
col.treatment	[character] column in the dataset indicating the treatment variable, that should be permuted within-individuals.
col.strata	[character] column in the dataset indicating the factor variable, to stratify the comparison on.
col.value	[character] column in the dataset indicating the value of the measurement.
col.cluster	[character] column in the dataset indicating the cluster level (e.g. patient identity).
n.perm	[numeric] number of permutation to be performed.
method.adj	[character] either "single-step" or "step-down" max adjustment. Both order the hypothesis to test based on the evidence shown by the data but the latter provides extra statistical power for rejecting the 2nd to last hypotheses.
seed	[integer,>0] set initial state of the random number generation (passed to set.seed).
trace	[logical] should a progress bar be display displaying the execution of the permutations.

Examples

library(mvtnorm)
library(reshape2)

## generate data
n <- 25

Sigma1 <- (diag(1-0.5,5,5)+0.5) * tcrossprod(c(1,2,1,2,3))
Sigma2 <- (diag(1-0.5,5,5)+0.5) * tcrossprod(c(1,2,1,2,3))
Sigma3 <- (diag(1-0.5,5,5)+0.5) * tcrossprod(c(1,2,1,2,3))
Sigma <- as.matrix(Matrix::bdiag(Sigma1,Sigma2,Sigma3))+0.25
mu1 <- c(1,2,1,2,1.5)
mu <- c(mu1, mu1+1.5,mu1-1)

set.seed(10)
dfW <- data.frame(1:n,mvtnorm::rmvnorm(n, mean = mu, sigma = Sigma))
colnames(dfW) <- c("id",paste0("score",1:5,"_B"),
                   paste0("score",1:5,"_P"),
                   paste0("score",1:5,"_K"))
dfL <- melt(dfW, id.vars = "id")
dfL$score <- sapply(strsplit(x=as.character(dfL$variable),split="_",fixed = TRUE),"[",1)
dfL$treatment <- sapply(strsplit(x=as.character(dfL$variable),split="_",fixed = TRUE),"[",2)

## test mean
t.test(dfW$score1_P,dfW$score1_B,paired = TRUE)
t.test(dfW$score5_P,dfW$score5_B,paired = TRUE)

if(require(MKinfer)){
perm.t.test(dfW$score1_P,dfW$score1_B,paired = TRUE)
perm.t.test(dfW$score2_P,dfW$score2_B,paired = TRUE)
perm.t.test(dfW$score3_P,dfW$score3_B,paired = TRUE)
perm.t.test(dfW$score4_P,dfW$score4_B,paired = TRUE)
perm.t.test(dfW$score5_P,dfW$score5_B,paired = TRUE)
}

resPerm <- pairedCompMean(dfL[dfL$treatment %in% c("P","B"),], n.perm = 1e3,
               col.treatment = "treatment", col.strata = "score", col.value = "value", col.cluster = "id",
               seed = NULL)
resPerm

bozenne/butils documentation built on Oct. 14, 2023, 6:19 a.m.