4
$\begingroup$

I regret that I cannot be overly explicit due to confidentiality agreements. I will try to state my problem by analogy.

Suppose that you have two groups of people, one with a condition and one normal. You make roughly 80 measurements of circumference of different body parts on each person. Clearly, these measurements will all be correlated, some more highly than others (measurements that are close to each other will tend to be more correlated than others).

You have an a priori hypothesis that the largest differences between groups will occur at certain places (e.g. "patients with the condition will have larger arms than patients without").

As a first step, we looked at t-tests comparing the measurements between the groups; there were quite a few significant differences and most of the largest differences were in the hypothesized region.

My idea was to use some form of permutation test to account for multiple testing while accounting for the correlations. This would not only be to test the specific hypothesis but to correct for multiple testing on all 80 measurements (instead of using Bonferroni). However, I need help:

1) Is this reasonable? 2) If reasonable, any guidance on how to do it in either R or SAS? (any hints or references welcome)

ADDING A REFERENCE One paper that does something similar to what I want is

Nichols and Holmes

But that is for FMRI data, which is somewhat different; and they give a program in Matlab, which I don't have or know. Also, they don't consider exactly the sort of problem I have.

$\endgroup$
1
+50
$\begingroup$

measurements that are close to each other will tend to be more correlated than others

If you have some kind of adjacency you could try looking for "clumps" (groups of adjacent interesting results) then comparing to what you get by either permuting the locations or group membership. It wouldn't necessarily have to be physical adjacency. Rather than choosing an arbitrary cutoff you can try to figure out the best case scenario in your favor and see how strong that is. Your statistic could be the size of the clump, the cutoff used, the percent of permutations of greater size at given cutoffs, or some index that is a combination. I've never seen this done anywhere but I am working on a similar problem.

The script below compares the largest clump in the image lower than a sequence of cutoffs in the range of .2 to .01 to those found by permutation results, then chooses the cutoff and clump with the lowest number of permutation results that were larger.

enter image description here

require(raster)
require(igraph)


image.plot(1:nrow(m), 1:ncol(m), m, breaks=seq(0,1,by=.01), col=rainbow(100), zlim=c(0,1))

alphas<-seq(.2,0.01,by=-.01)
min.size<-3


real.data=matrix(nrow=length(alphas),ncol=2)
for(i in 1:length(alphas)){
  temp.m<-m
  temp.m[which(temp.m>alphas[i])]<-NA
  if(length(which(is.finite(temp.m)))<min.size){
    real.data[i,]<-cbind(alpha,0)
  }else{
    clumps<-clump(raster(temp.m),directions=8)
    clump.data<-freq(clumps, useNA="no")
    biggest.clump<-max(clump.data[,2])
    real.data[i,]<-cbind(alphas[i],biggest.clump)
  }
}
colnames(real.data)<-c("alpha","Clump.Size")



permutations<-200
sim.data=matrix(nrow=permutations,ncol=length(m))
for(i in 1:permutations){
  sim.data[i,]<-sample(m, length(m), replace=F)
}


sim.clump=matrix(nrow=length(alphas)*nrow(sim.data),ncol=2)
cnt<-1
for(i in 1:length(alphas)){
  for(p in 1:nrow(sim.data)){
    temp.m<-m
    temp.m[1:length(temp.m)]<-sim.data[p,]
    temp.m[which(temp.m>alphas[i])]<-NA
    if(length(which(is.finite(temp.m)))<min.size){
      sim.clump[cnt,]<-cbind(alphas[i],0)
      cnt<-cnt+1
    }else{
      clumps<-clump(raster(temp.m),directions=8)
      clump.data<-freq(clumps,useNA="no")
      biggest.clump<-max(clump.data[,2])
      sim.clump[cnt,]<-cbind(alphas[i],biggest.clump)
      cnt<-cnt+1
    }
  }
}

compare=matrix(nrow=nrow(real.data),ncol=3)
for(i in 1:nrow(real.data)){
  num<-length(which(sim.clump[,1]==real.data[i,1] & sim.clump[,2]>=real.data[i,2]))
  percent<-num/nrow(sim.data)
  compare[i,]<-c(real.data[i,],percent)
}
colnames(compare)<-c("alpha","Top.Clump.Size","Percent")

best.clump<-compare[which(compare[,3]==min(compare[,3])),]

if(length(best.clump)>3){
  best.clump<-best.clump[which(best.clump==min(best.clump[,1])),]
}


temp.m<-m
temp.m[which(temp.m>best.clump[1])]<-NA
clump.map<-as.matrix(clump(raster(temp.m)))
best.clump.id<-as.numeric(names(table(clump.map)[1]))
best.clump.coords<-which(clump.map==best.clump.id,arr.ind=T)


dev.new()
image.plot(1:nrow(m),1:ncol(m), m, breaks=seq(0,1,by=.01), col=rainbow(100),    
zlim=c(0,1))
points(best.clump.coords,pch=16)


> dput(m)
structure(c(0.865102812502504, 0.836453341655039, 0.90370356257473, 
0.920150007969666, 0.834252821321425, 0.88691326798952, 0.95625021357218, 
0.966077554549026, 0.876618617922488, 0.800510067385609, 0.140714584096359, 
0.165367840985583, 0.887065809921221, 0.88207475582995, 0.869682335922427, 
0.181837532928766, 0.815218422415644, 0.90035552309918, 0.890630877655364, 
0.162942811596665, 0.162475967113964, 0.814938152638217, 0.189192052732397, 
0.808682406205577, 0.184229286652596, 0.166204898316331, 0.164311765529674, 
0.852210918877963, 0.821143172068156, 0.871165208495517, 0.159637928956964, 
0.821378863048212, 0.171458731465239, 0.109475889431187, 0.0847341126305428, 
0.150613887165366, 0.170171167245293, 0.17399055067852, 0.813975277166275, 
0.840917132279195, 0.871500443905951, 0.863732557825426, 0.127743512687412, 
0.0841117221895071, 0.0856239370342532, 0.165967905516149, 0.1664115976353, 
0.82218002302312, 0.846348744336673, 0.877288700874692, 0.952919232078247, 
0.853705806889241, 0.181541385029858, 0.127547084111945, 0.170999488382186, 
0.154020138365571, 0.135911184927553, 0.186352521809119, 0.839313542020112, 
0.858507627504672, 0.965494682207016, 0.816722462771706, 0.174153259912455, 
0.181018743063141, 0.165433887212926, 0.195108282530258, 0.80659485948176, 
0.09554158394495, 0.106164217572178, 0.0799071434463993, 0.930029405018248, 
0.820277463053794, 0.151569946447117, 0.131110764891388, 0.177503429054507, 
0.83846863096766, 0.17359739352256, 0.0602178345535283, 0.0668336042151074, 
0.0491473775305895, 0.948812369805206, 0.814305631855217, 0.123005923377645, 
0.139738001159334, 0.831513371125685, 0.198732850276339, 0.117688539783991, 
0.0607235595181853, 0.0465879440856685, 0.185091314291344, 0.863157732274661, 
0.118637075448383, 0.0916003666812711, 0.84342604027262, 0.889141857102965, 
0.865196613575319, 0.819165829895222, 0.118277219263935, 0.130642430905136, 
0.139540111909131), .Dim = c(10L, 10L))

enter image description here

> best.clump
         alpha Top.Clump.Size        Percent 
          0.07           5.00           0.00 
$\endgroup$
2
  • $\begingroup$ That is very cool. I am going to give it the reward. $\endgroup$ – Peter Flom Nov 27 '13 at 10:52
  • $\begingroup$ @PeterFlom I was thinking that this method could also incorporate prior probability by creating a map of weights/ prior probabilities that adjusts the values in the map of the data. $\endgroup$ – Flask Nov 27 '13 at 15:42
1
$\begingroup$

The first issue that appears to me is that of multiplicity. If our objective is to consider each of the 80 places on the body where these individuals are measured and test whether there's a difference between diseased and healthy individuals, then having a 0.05 level test at each of the 80 places will result in an average of 4 type I errors. That's no good.

Correcting this by using Bonferroni would yield a valid test for each site but is extremely conservative due to the assumed independence between each anatomical site (which we know to be false).

If $n \gg 80$, then it would be possible to simultaneously test the partial correlations for conditional independence using a multivariate regression model with fixed effects for each site and an interaction with disease indicator. Using an L-1 penalty for the interaction parameters would force small differences toward zero and give you more power to estimate larger differences. Alternately, if data were matched (paired), you can compute differences in diseased and undiseased individuals and fit the differences model with intercept through the origin.

If you're doing a randomization / permutation, it should be possible to calibrate it using a simulation. My question to you is this: the permutation test will only give you the joint sampling distribution of the model parameters under the null hypothesis. What next? What's your plan to use that information and obtain a calibrated and well powered test of hypothesis about the 80 sites?

$\endgroup$
4
  • $\begingroup$ First, N is not >> 80. It's about 120. However, it is matched on age and sex. What differences model are you suggesting? I would want the permutation test to simply say "The chance that these differences would occur, given these correlations but no difference between control and treatment is XXX" I know it will be very small (having looked at the data) but I need something more than "Ecce! Differences!" $\endgroup$ – Peter Flom Nov 19 '13 at 23:23
  • $\begingroup$ Matching on age and sex is a start, but generally I take only data which are matched upon extremely prognostic factors to be really treated as matched data. (In this case, we would like to know perhaps BMI, body fat percentage, urinalysis and other factors causally related to outcome and potentially correlated with exposure). The differences model is described nicely by Carlin, 2005 and I only mention it since there were equal numbers of cases and controls. $\endgroup$ – AdamO Nov 20 '13 at 17:55
  • $\begingroup$ The problem that remains for me is exactly how the permutation test will be carried out to handle multiple testing. This is new to me. $\endgroup$ – AdamO Nov 20 '13 at 17:56
  • $\begingroup$ There has been work on this, using single threshold and supra threshold; I will post some in just a bit. $\endgroup$ – Peter Flom Nov 20 '13 at 17:58
1
$\begingroup$

I think your permutation idea is solid. Assume there are only two measurements y1 and y2. Here's how I'd do it in R:

set.seed(1245)
group <- c(rep(0, 50), rep(1, 50))
z <- rnorm(100)  # Direct common cause of y1 and y2
y1 <- .7 * group + 1.5 * z + rnorm(100) # explanitory variable
y2 <-  0 * group + 1.5 * z + rnorm(100)

t1.stat <- t.test(y1[group == 0], y1[group == 1])$statistic
    t2.stat <- t.test(y2[group == 0], y2[group == 1])$statistic

cat("t statistics - Var 1:", t1.stat, " Var 2:", t2.stat, "\n")

t1.collection <- c()
t2.collection <- c()

for( rep in 1:1000) {
  group <- sample(group, 100)
  t1 <- t.test(y1[group == 0], y1[group == 1])$statistic
  t2 <- t.test(y2[group == 0], y2[group == 1])$statistic
  t1.collection <- c(t1.collection, t1)
  t2.collection <- c(t2.collection, t2)
}
hist(t1.collection)
hist(t2.collection)

cat("Permutation p-values - Var 1 ", mean(t1.collection >= abs(t1.stat)),
    "Permutation p-values - Var 2 ", mean(t2.collection >= abs(t2.stat)), "\n")

The one thing that this doesn't cover is that you have prior beliefs about which variables are most likely to pop. Sounds like you'd need to go bayesian for that, and your model looks like a standard multivariate regression (all variables have the same single group regressor). Sometimes I wish I was a bayesian!

$\endgroup$
3
  • $\begingroup$ Thanks for this; but, as I said in the post, there are about 80 measurements and they are all correlated to different degrees. Also, in the definition of y1 and y2, where did 0.7 and 0.0 come from? $\endgroup$ – Peter Flom Nov 22 '13 at 18:05
  • $\begingroup$ Yeah, you're right - I didn't fully answer the question, at least not the multiple testing component. I thought that two variables would be a good place to start, though you have 80 in actuality. 0.7 is arbitrary and signifies a true discovery. The 0.0 is meant to make it explicit that there is no group effect for y2 and it would be a false discovery if found. Thinking about the permutation based method, it always keeps the rows together, even as it permutes the group labels - I guess I thought that it was "preserving the correlation structure," but perhaps I didn't think it through. $\endgroup$ – Ben Ogorek Nov 22 '13 at 23:21
  • $\begingroup$ Thanks! I did upvote your answer, but not the +50. $\endgroup$ – Peter Flom Nov 22 '13 at 23:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.