Pairs analyses? #54
Description
Hi @emhart @ngotelli, hope you both had nice weekends.
My Windows box doesn't want to run Dr. Ulrich's binary (it thinks the .exe file might be a virus for some reason). I was wondering if either of you would be interested in looking at my implementation (below) and letting me know if it looks correct to you. If it does, I'd be happy to submit it as a pull request.
Thanks!
library("EcoSimR")
set.seed(1)
# create “observed” data matrix ----------------------------------------------
m = matrix(rbinom(50^2, size = 1, prob = .1), ncol = 100)
n_samples = 999 # number of samples to collect
thin = 1000 # number of iterations between samples
n_species = nrow(m)
# I think this is a correct way to get scaled pairwise C-scores...
scaled_pairwise_c_score = function(x){
shared = tcrossprod(m)
sums = rowSums(m)
upper = upper.tri(shared)
raw_scores = (sums[row(shared)[upper]] - shared[upper]) *
(sums[col(shared)[upper]] - shared[upper])
raw_scores / tcrossprod(sums)[upper]
}
# Record observed c_scores
observed = scaled_pairwise_c_score(m)
# Null --------------------------------------------------------------------
# Create empty matrix to store null results for all n-choose-2 pairs
null = matrix(NA, nrow = choose(n_species, 2), ncol = n_samples)
for(i in 1:n_samples){
for(j in 1:thin){
# Update a
m = sim9_single(m)
}
# Fill in the next column of the null distribution
null[ , i] = scaled_pairwise_c_score(m)
}
# Calculate P-values ------------------------------------------------------
# For each observed value, compare it to the null distribution from _ALL_
# species pairs to calculate p-values
p = sapply(observed, function(x) mean(x < null))