# Probability of multiple parasites given the presence of one parasite

I am trying to simulate parasitic plants that infect trees. I would like to test wether the presence of one parasite makes it more likely that the tree has more than one parasite.

The number of trees with n parasites is

10822   887   375   168    74    25    12     6     4     1


meaning that 10822 trees have 0 parasites, 887 trees have 1 parasite and so forth.

I assumed that if the parasites would be distributed randomly then they would be described by a Poisson distribution. First question is: is this reasonable?

If it is then I should test this hypothesis by comparing the distribution with a Poissonian.

This is what I did in R

set.seed(5)
y= c(10822,887,375,168,74,25,12,6,4,1)
# Calculate expected value
lambda = weighted.mean(x=seq(0,length(y)-1), w=y)
# Generate Poisson distribution with mean lambda
x.poi = table(rpois(n=sum(y), lambda = lambda))
# Calculate expected frequencies
ex.freq = as.vector(x.poi/sum(x.poi))
# Only frequencies with at least 5% prob. Group the others
while(ex.freq[length(ex.freq)] < 0.05){
ex.freq[length(ex.freq) - 1] = ex.freq[length(ex.freq) - 1] + ex.freq[length(ex.freq)]
ex.freq = ex.freq[-length(ex.freq)]
}
# Club all observed outcomes outside of the range of expected frequencies so in my case since ex.freq is of size 3
# y -> c(10822, 887, 665)
y[length(ex.freq)] = sum(tail(y,length(y) - length(ex.freq) + 1))
# Remove other elemets
y = y[-seq(length(ex.freq) + 1,length(y))]
# Perform Chi Squared Test
chisq.test(y,p=ex.freq)


What I get in this case is

    Chi-squared test for given probabilities

data:  y
X-squared = 1348.2, df = 2, p-value < 0.00000000000000022


I am not sure how to interpret this data. Does this mean that my data is not distributed according to a Poisson? Or that it is?

EDIT

So as some of you asked, I'll try to give some detail about the processes that I am simulating. There are a few different scenarios so I'll try to be as clear as I can. The parasites I talk about are woody vines (lianas). Upon recruitment they are assigned to hosts in two ways:

1) Distributed at random only to trees.

2) Distributed at random among both trees and vines that are already attached to a tree.

Given this mode of recruitment I would expect the distribution to be Poissonian in the first case and non-Poissonian in the second case. I tested this with the above code and it seems to be confirmed.

Now if I let my forest grow, the initial distribution will change.

If I make a test at the quasi-equilibrium state, it seems that with both recruitment modes 1 and 2, the final distribution is non-Poissonian. So that would mean that during growth, there arre effects that remove the randomness and select for a more clumped distribution. No process directly affects the vines, so this has to be a tree-mediated effect. The more vines a tree has, the higher mortality it will experience.

• If your chi-squared test is correct, the very small P-value indicates rejection of the null hypothesis that the data are Poisson with the estimated $\lambda.$ However, I am not sure that the chi-squared test was done correctly. // Also, you have many more 0-counts than seem consistent with Poisson. Perhaps some trees in the forest are outside of the area of parasite infestation. – BruceET Sep 13 '19 at 18:03
• My idea was to actually prove that trees with one parasite are more likely to have more than one parasite. So in fact if the data are not Poisson distributed this would confirm my hypothesis. Is this line of reasoning correct? – Manfredo Sep 14 '19 at 14:39
• Why do you think the chi-squared was not done correctly? Is there something I should change? – Manfredo Sep 14 '19 at 14:43
• (a) There are many non-Poisson distributions, so it is hard to conclude what you want just from non-Poisson. Your conjecture is reasonable, but I don't see how this establishes it. (b) I do not follow the logic of your (uncommented) R code. How about a brief discussion in English about the procedure, which is unfamiliar --at least to me. In particular, where is your estimated $\lambda?$ why are you generating data? specifically, what are your expected frequencies under $H_0?$ why do you think c(10822, 887, 665) fairly represents your data? – BruceET Sep 14 '19 at 19:13
• Hi Manfredo, I like your problem and I have an idea on how to attack it. Is this for a scientific reaserch? Let me know. – Nicola Mingotti Sep 15 '19 at 22:00

First, your data is probably not distributed like Poisson. One of the first things we learn about the Poisson distribution is that mean equals variance. So the first to do with your data is:

lambda <- weighted.mean(0:9, y)
varhat <- Hmisc::wtd.var(0:9, y)
c(lambda, varhat)
[1] 0.2195733 0.4990236


so the estimated variance is more than double the mean. This indicated your data is overdispersed, rather common with count data (search this site for many examples.) Continuing, look at the estimated probabilities assuming a Poisson distribution, compared to the observed relative frequencies:

phat <- y/sum(y); names(phat) <- 0:9
round(rbind(phat, dpois(0:9, lambda)), 2)
0    1    2    3    4 5 6 7 8 9
phat 0.87 0.07 0.03 0.01 0.01 0 0 0 0 0
0.80 0.18 0.02 0.00 0.00 0 0 0 0 0


so it is quite clear this is not a Poisson distribution (given your large sample.) This should not be a surprise, overdispersion is a very common phenomenon! also in biostatistics.

Your data is about parasitic plants that infect trees. Do you know anything about how this parasitic plants spread? That should be your point of departure now. By the roots ... blowing in the wind ... tell us about that!