1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677 |
- function [obsDiff, pVal] = nestetPermutationTest(sample1, sample2, ...
- groupsCell, reorderedFlyIndsCell, ...
- permutations, sidedness, plotresult)
- allobservations = [sample1, sample2];
- observeddifference = nanmean(groupsCell{1}) - ...
- nanmean(groupsCell{2});
- effectsize = observeddifference / ...
- nanmean([nanstd(groupsCell{1}), nanstd(groupsCell{2})]);
- exact = 0;
- showprogress = permutations / 10;
- for n = 1:permutations
- if showprogress && mod(n,showprogress) == 0, display(sprintf('Permutation %d of %d', n, permutations)); end
- % selecting either next combination, or random permutation
- if exact, permutation = [allcombinations(n,:), setdiff(1:numel(allobservations), allcombinations(n,:))];
- else, permutation = randperm(length(allobservations)); end
- % dividing into two samples
- randomSample1 = cell(1, numel(sample1));
- randomSample2 = cell(1, numel(sample2));
- jCount = 1;
- for iElem = permutation
- % First fill sample1, then check if the sample coes from first or
- % second genotype.
- if jCount <= length(sample1)
- if iElem > length(sample1)
- sampleInds = reorderedFlyIndsCell{2} == iElem - length(sample1);
- randomSample1{jCount} = groupsCell{2}(sampleInds);
- else
- sampleInds = reorderedFlyIndsCell{1} == iElem;
- randomSample1{jCount} = groupsCell{1}(sampleInds);
- end
- else
- if iElem > length(sample1)
- sampleInds = reorderedFlyIndsCell{2} == iElem - length(sample1);
- randomSample2{jCount - numel(sample1)} = groupsCell{2}(sampleInds);
- else
- sampleInds = reorderedFlyIndsCell{1} == iElem;
- randomSample2{jCount - numel(sample1)} = groupsCell{1}(sampleInds);
- end
- end
- jCount = jCount + 1;
- end
- % saving differences between the two samples
- randomdifferences(n) = nanmean(cat(2, randomSample1{:})) - ...
- nanmean(cat(2, randomSample2{:}));
- end
- % getting probability of finding observed difference from random permutations
- if strcmp(sidedness, 'both')
- p = (length(find(abs(randomdifferences) > abs(observeddifference)))+1) / (permutations+1);
- elseif strcmp(sidedness, 'smaller')
- p = (length(find(randomdifferences < observeddifference))+1) / (permutations+1);
- elseif strcmp(sidedness, 'larger')
- p = (length(find(randomdifferences > observeddifference))+1) / (permutations+1);
- end
- % plotting result
- if plotresult
- figure;
- hist(randomdifferences, 20);
- hold on;
- xlabel('Random differences');
- ylabel('Count')
- od = plot(observeddifference, 0, '*r', 'DisplayName', sprintf('Observed difference.\nEffect size: %.2f,\np = %f', effectsize, p));
- legend(od);
- end
- obsDiff = observeddifference;
- pVal = p;
- end
|