Inequalities for tail extreme tail probabilities not correct with RA1

[ngotelli]

Need to check code for summary function in niche overlap with ra1and make sure it is handling extreme cases the same way for all algorithms:

# If observed index is smaller than any of the 1000 simulated values, 
# tail inequalities should read:
p(Obs <= Null) < 0.001
p(Obs >= Null) > 0.999

[emhart]

@ngotelli Is this not correct?

> warbMod <- niche_null_model(dataMacWarb,algo="ra1",metric="pianka",nReps=10000)
  |=============================================================================================| 100%
> summary(warbMod)
Time Stamp:  Thu Feb 26 22:26:14 2015 
Reproducible:  
Number of Replications:  
Elapsed Time:  2.5 secs 
Metric:  pianka 
Algorithm:  ra1 
Observed Index:  0.55514 
Mean Of Simulated Index:  0.75565 
Variance Of Simulated Index:  0.001333 
Lower 95% (1-tail):  0.69415 
Upper 95% (1-tail):  0.81361 
Lower 95% (2-tail):  0.68064 
Upper 95% (2-tail):  0.82395 
P(Obs <= null) =  0 
P(Obs >= null) =  1 
P(Obs = null) =  0 
Standardized Effect Size (SES):  -5.4919 
> plot(warbMod)
> min(warbMod$Sim)
[1] 0.5928378
> warbMod$Obs
[1] 0.5551383

The p-values look right in this example.

[ngotelli]

Hi @emhart

This is all correct except the printout has to be modified in this case because the observed is more extreme than any of the simulated values. With 10,000 replicates (oddly, the 10,000 is not displayed in this summary after Number of Replications: ) the output should be:

P(Obs <= null) < 0.0001
P(Obs >= null) > 0.9999
P(Obs = null) = 0

and if the analysis was based on only 10 replications, it should print:

P(Obs <= null) < 0.1
P(Obs >= null) > 0.9
P(Obs = null) = 0

When the observed index value is within the range of the simulated values, then the simple frequentist calculation is fine. But once it is outside, we want to display it this way to eliminate probability estimates of 1.0 or 0.0, which are biased because those values are not possible for a true stochastic event. It also means that the more extreme the actual probability, the more simulations you need to estimate it properly. Otherwise, the inference is just a region of probability, not a point estimate. So the summary code has this tripartite fork to handle the two extreme cases and the more typical intermediate case.

I had intended to clean all of these up, and I think I may already have done so for co-occurrence.

thanks,

N.

[emhart]

@ngotelli I was looking more closely at the code for your if statements, and I just want to go over the logic. There's three different conditionals in the summary statement. I pulled them all out into a summary function that evaluates each and then shows the output of EcoSimR. I've labeled them case 1 2 and 3 (for the order they are in the if statement). The main issues are 1). case 1 and 2 are the same conditional, but because it's an if else statement case 2 is never actually evaluated. I lay out three scenarios, obs >>>> simulation, obs ~= sim, and obs <<<< sim. I think scenario 1 gives the wrong output. Can you check that my thinking is correct and not subject to late night parent brain? Thanks.

nullEval <- function(sim,obs){

 case1 <- obs > max(sim)
 case2 <- sum(obs <= sim)==0
 case3 <- (case1 + case2) == 0
 case1.pval <- c(1/length(sim),(length(sim) - 1)/length(sim))
 case2.pval <- c((length(sim) - 1)/length(sim),1/length(sim))
 case3.pval <- c(sum(obs >= sim)/length(sim),sum(obs <= sim)/length(sim))

 cat("Case 1 is ",case1,"\n")
 cat("Case 1 P(Obs <= null)", case1.pval[1],"\n")
 cat("Case 1 P(Obs >= null)", case1.pval[2],"\n")

 cat("Case 2 is ", case2,"\n")
 cat("Case 2 P(Obs <= null)", case2.pval[1],"\n")
 cat("Case 2 P(Obs >= null)", case2.pval[2],"\n")
 cat("Case 3 is ", case3,"\n")
 cat("Case 3 P(Obs <= null)", case3.pval[1],"\n")
 cat("Case 3 P(Obs >= null)", case3.pval[2],"\n")

  cat("######################","\n")
 cat("What EcoSimR would output","\n")
  if (case1) {
    cat("P(Obs <= null) < ",1/length(sim),  "\n")
    cat("P(Obs >= null) > ",(length(sim) - 1)/length(sim),  "\n")
  } else if (case2) {
    cat("P(Obs <= null) > ",(length(sim) - 1)/length(sim), "\n")
    cat("P(Obs >= null) < ", 1/length(sim), "\n")
  } else {
    cat("P(Obs <= null) = ", format(sum(obs >= sim)/length(sim),digits=5),  "\n")
    cat("P(Obs >= null) = ", format(sum(obs <= sim)/length(sim),digits=5), "\n")
  }

}

Scenario 1

Observed value > than all simulated values In this case, I would expect that probability of observed > null to be ~1 and P obs < null ~0

sim <- rnorm(10000,10,2)
obs <- 100
nullEval(sim,obs)
Case 1 is  TRUE 
Case 1 P(Obs <= null) 1e-04 
Case 1 P(Obs >= null) 0.9999 
Case 2 is  TRUE 
Case 2 P(Obs <= null) 0.9999 
Case 2 P(Obs >= null) 1e-04 
Case 3 is  FALSE 
Case 3 P(Obs <= null) 1 
Case 3 P(Obs >= null) 0 
###################### 
What EcoSimR would output 
P(Obs <= null) <  1e-04 
P(Obs >= null) >  0.9999 

Note that both conditionals in case 1 and 2 evaluate to true, but the outputs are opposite. Case 1 is printed, however wouldn't we expect case 2 to be the correct output?

Scenario 2

Observed is close to null mean

sim <- rnorm(10000,10,2)
obs <- 10
nullEval(sim,obs)

Case 1 is  FALSE 
Case 1 P(Obs <= null) 1e-04 
Case 1 P(Obs >= null) 0.9999 
Case 2 is  FALSE 
Case 2 P(Obs <= null) 0.9999 
Case 2 P(Obs >= null) 1e-04 
Case 3 is  TRUE 
Case 3 P(Obs <= null) 0.5005 
Case 3 P(Obs >= null) 0.4995 
###################### 
What EcoSimR would output 
P(Obs <= null) =  0.5005 
P(Obs >= null) =  0.4995 

Here case 3 evaluates to true, and both 1 and 2 are false. The p-values are ~ 0.5 each, which is what I would exect.

Scenario 3

Observed value < than all simulated values In this case, I would expect that probability of observed > null to be ~0 and obs < null ~1

Here's what we output:

  sim <- rnorm(10000,10,2)
  obs <- -100
 nullEval(sim,obs)

Case 1 is  FALSE 
Case 1 P(Obs <= null) 1e-04 
Case 1 P(Obs >= null) 0.9999 
Case 2 is  FALSE 
Case 2 P(Obs <= null) 0.9999 
Case 2 P(Obs >= null) 1e-04 
Case 3 is  TRUE 
Case 3 P(Obs <= null) 0 
Case 3 P(Obs >= null) 1 
###################### 
What EcoSimR would output 
P(Obs <= null) =  0 
P(Obs >= null) =  1 

Again case 3 is true, and the output is what I would expect.

So it seems that case 1 is never right. It's also strange that conditionals in case 1 and 2 are just different ways of writing the same conditional, but you have the outputs opposite in each.

[ngotelli]

Hi @emhart:

Thanks for checking this and setting up this useful nullEval function. Yes, there was a problem in the fork for Case2, but there was also a problem in the calculation of the tail probabilities. In brief, to get the tail probability for the inequalities, you actually need to get the count of extreme values in the other tail.

When it is working properly, obs values that are very large compared to sim will trigger Case1. This will give you an upper tail p very close to 0, and lower tail p very close to 1, obs > mean(sim), and a SES (standardized effect size = (obs - mean(sim)/sd(sim)) with a large positive value.

obs values that are very small compared to sim will trigger Case2. This will give you an upper tail p very close to 1, and lower tail p very close to 0, obs < mean(sim), and a SES (standardized effect size = (obs - mean(sim)/sd(sim)) with a large negative value.

Values that are not extreme will trigger Case3. If obs is very close to mean(sim), the upper and lower p values will both be close to 0.5, and the SES will be close to 0. The larger obs is than mean(sim), the more positive the SES, and the smaller obs is than mean(sim), the more negative the SES.

So, here are all the corrections to the nullEval function. I added obs, mean(sim), and SES so you could quickly make these comparisons.

# Corrected nullEval function

nullEval <- function(sim,obs){

 case1 <- obs > max(sim)
 case2 <- obs < min(sim)
 case3 <- (case1 + case2) == 0
 case1.pval <- c((length(sim) - 1)/length(sim),1/length(sim))
 case2.pval <- c(1/length(sim),(length(sim) - 1)/length(sim))
 case3.pval <- c(sum(obs >= sim)/length(sim),sum(obs <= sim)/length(sim))

 cat("Case 1 is ",case1,"\n")
 cat("Case 1 P(Obs <= null)", case1.pval[1],"\n")
 cat("Case 1 P(Obs >= null)", case1.pval[2],"\n")

 cat("Case 2 is ", case2,"\n")
 cat("Case 2 P(Obs <= null)", case2.pval[1],"\n")
 cat("Case 2 P(Obs >= null)", case2.pval[2],"\n")
 cat("Case 3 is ", case3,"\n")
 cat("Case 3 P(Obs <= null)", case3.pval[1],"\n")
 cat("Case 3 P(Obs >= null)", case3.pval[2],"\n")

 cat("obs = ", format(obs),"\n")
 cat("mean(sim)", format(mean(sim),digits=5),"\n")
 cat("SES = ", format((obs - mean(sim))/sd(sim),digits=5),"\n")

  cat("######################","\n")
 cat("What EcoSimR would output","\n")
  if (case1) {
    cat("P(Obs <= null) > ",(length(sim) - 1)/length(sim),  "\n")
    cat("P(Obs >= null) < ",1/length(sim),  "\n")
  } else if (case2) {
    cat("P(Obs <= null) < ", 1/length(sim), "\n")
    cat("P(Obs >= null) > ",(length(sim) - 1)/length(sim), "\n")
  } else {
    cat("P(Obs <= null) = ", format(sum(obs >= sim)/length(sim),digits=5),  "\n")
    cat("P(Obs >= null) = ", format(sum(obs <= sim)/length(sim),digits=5), "\n")
  }
}

Run a few test cases and see if you agree. If this looks OK, I'd be grateful if you could drop in the corrected lines code into the correct places for all of the EcoSimR scripts that use it.

Thanks, and sorry for the delay in getting back to you on this one.

[emhart]

Looks good to me, I'll integrate the code tonight. I noticed however that you swapped the probabilities calculation around. In the old version it was:

 cat("P(Obs <= null) < ",1/length(sim),  "\n")
 cat("P(Obs >= null) > ",(length(sim) - 1)/length(sim),  "\n")

In the new version you put

    cat("P(Obs <= null) > ",(length(sim) - 1)/length(sim),  "\n")
    cat("P(Obs >= null) < ",1/length(sim),  "\n")

Why did we end up switching it? I ask because my intuition is that if the obs = 1 and the null = 0, then the P obs >= null ~ 1 (e.g. the probability that 1 is greater than a null with a mean of 0) However P null>=obs ~0 or 1 - (p obs>=null).

If the question is P(obs == null), then I'd say that's also 1 - (p obs>=null). Am I incorrect in my thinking?

Perhaps it's just an issue of semantics. but if Obs >>> Null distribution, how can the probability that the Obs is > Null be 0?

[ngotelli]

This is definitely an issue of confusing semantics! For conciseness in the output, we print P(obs > sim). But what we are actually calculating is P(a randomly drawn value from the null distribution is as large or larger than the observed value). As you note, if Obs >>>> Null then it seems nearly certain that P(obs > sim) ~ 1.0 because most of the probability mass is to the left of Obs. But if Obs was actually sampled from Null, the P(Obs >>>>> Null) is estimating the size of the upper tail of the distribution. That tail is very small when Obs is large. I realize that it might be more logical to describe this as P(null > obs), but the short-hand P(Obs > Null) is following long-standing convention on how P-values from simulations are reported in the literature.

If you look at the hist output plot, you will see that the area under the curve above and below the red line (= obs) should match the probability calculations of the upper tail (P(Obs > Null)) and the lower tail (P(Obs < Null)) respectively, even though they seem opposite of the short-hand notation.

Hope that helps and doesn't make it even more confusing!

I believe the new code is correctly calculating upper- and lower-tail probabilities for all possible ranges of Obs. It is just a question of what we call them. Aaron, any thoughts on this?

[amellison17]

On the one hand, we should stick with convention in the literature. But what we're really looking at is the tail probability P(Obs < Null). Originally, the assumption wasthat competition structures assemblages, and so Obs should always be way out to the left. It didn't make sense in the 70s and 80s to ever think that Obs would ever be greater than null. Now we know better. But in the context of the underlying hypothesis, I always think of it as :P(Obs < Null).

Why don't we just write it as P(Obs)?

[ngotelli]

I don't think we can say P(Obs) because we are reporting both tails: P(Obs <= null) and P(Obs >= null). If Obs never matches any of the values in null, then of course these sum to 1.0. We need to report both tails because, depending on the metric and the question, either or both tails may be of interest. Also, we don't want to say P(Obs), because we already report that, which is the probability of an exact match.

[amellison17]

Then it probably should be P (Obs < Null) because that is what most people are interested in, I think.

[ngotelli]

But it still depends on the metric! For species segregation, that would be indicated by a high C-score in co-occurrence (upper tail), but a low pianka score in niche overlap (lower tail). What if we reported:

Lower-tail P = 0.971
Upper-tail P = 0.029

and for the extreme cases

Lower-tail P > 0.999
Upper-tail P < 0.001

The only thing I don't like about this is we are not reporting the ties. Perhaps we need one additional output, just in terms of the raw frequencies of the different cases:

Observed metric > 971 simulated metrics
Observed metric <  29 simulated metrics
Observed metric = 0 simulated metrics

If we report :"Upper-tail p", "Lower-tail p" and these 3 frequencies, does that make it crystal clear?

[amellison17]

that's much better.

[ngotelli]

OK. Ted, can you please set it up that way for the cat statements in all of the summary functions? Thanks.