# Are these populations non-random and different?

What would be the simplest, straightest-forwardest way to determine the following:

1. Whether A,B,C is non-randomly distributed in each single group (IE: are there more A's in group 1 than random chance would predict?)
2. Whether the distribution of A,B,C is different between the groups. (IE: Group 1 is has more Cs than Bs, so does Group 2 - but is the difference between the two groups significant?)

Data:

Group,  A,   B,   C,
1,   126, 357, 348   N=382
2,   86,  196, 139   N=207
3,   63,  185, 162   N=193
4...
5...


A <- c(126,86,63,54,47,40,32,32,29,29,27,26,20,18,14)
B <- c(357,196,185,137,95,74,45,69,64,49,54,80,62,41,56)
C <- c(348,139,162,126,82,69,35,63,40,42,40,55,44,29,35)
N <- c(382,207,193,143,100,80,45,70,70,53,55,84,67,42,57)


A,B,C are counts of presence/absence of traits within each group, so lots of As in group 1 are not . Groups are each separate non-related populations.

Example: lets say each group represents a species of lizard, while A, B, C indicate whether the lizard has spots on its head(A), Body(B), or Tail(C). For species 1 (Group 1) 382 lizards were examined, 126 had spots on their head, 357 on their body, 348 on their tail... Of 207 lizards in species 2, 86 had spotty heads, 196 spotty bodies, 139 spotty tails.

So, is spottiness non-random for members in a group? And does spottiness vary significantly between the groups?

I think this is a basic fundamental question, but while I have been bashing through countless pages of theory explaining different models looking at more complex situations, I have totally lost sight of the basics -- which I have only just recently learned.

Thanks for helping me get back on track.

-

For question 2, I would start by looking at the distributions of the percentages of each group that have a given trait. e.g.:

A <- c(126,86,63,54,47,40,32,32,29,29,27,26,20,18,14)
B <- c(357,196,185,137,95,74,45,69,64,49,54,80,62,41,56)
C <- c(348,139,162,126,82,69,35,63,40,42,40,55,44,29,35)
N <- c(382,207,193,143,100,80,45,70,70,53,55,84,67,42,57)

A <- A/N
B <- B/N
C <- C/N


Then you can make some stem and leaf plots to examine the distributions:

> stem(A)

The decimal point is 1 digit(s) to the left of the |

2 | 5
3 | 01338
4 | 123679
5 | 05
6 |
7 | 1

> stem(B)

The decimal point is 2 digit(s) to the left of the |

90 | 4
92 | 5555
94 | 70289
96 | 6
98 | 226
100 | 0

> stem(C)

The decimal point is 1 digit(s) to the left of the |

5 | 7
6 | 15679
7 | 389
8 | 2468
9 | 01


All three look somewhat normal, but that's pretty subjective. Only the .71 in stem(A) appears to be an outlier to me.

-
Nice way to visualize. Thanks. But how would you quantify this? Maybe there is not a way. Don't get me wrong, I'm not trying to be difficult, I just have no idea. –  Trees4theForest Dec 16 '11 at 4:08

I'm not sure I fully understand the situation (your data or your question). @Zach has some good ideas, so I thought I'd follow his lead and throw out some information and we'll see if something helps.

1) To determine if data are random, the runs test can be used. A run is a series of data points that are similar in some way. There are (at least) two ways of doing a runs test. First, you can check if each new data point is higher or lower than the previous one, alternatively, you can check if each data point is high or low relative to the median. Both of these assume that your data are ordered by the temporal sequence in which they were generated. From there, you use one of the methods above to convert your data into a string of ones and zeros. Sequential values of the same type are a run (e.g., 1101 is 3 runs). The number of runs is a binomial random variable. Here is some code for R:

# compute proportions & create vectors for runs
Aprop   = A/N;     Bprop   = B/N;     Cprop   = C/N
Achange = c();     Bchange = c();     Cchange = c()
Ahigh   = c();     Bhigh   = c();     Chigh   = c()

for(i in 1:length(N)) {
if(i>1) {
# determine if values went up or down
Achange[i-1] = (Aprop[i]-Aprop[i-1])>0
Bchange[i-1] = (Bprop[i]-Bprop[i-1])>0
Cchange[i-1] = (Cprop[i]-Cprop[i-1])>0
}
# determine if values are high or low
Ahigh[i] = Aprop[i]>median(Aprop)
Bhigh[i] = Bprop[i]>median(Bprop)
Chigh[i] = Cprop[i]>median(Cprop)
}

# conduct analyses
round(Aprop, 2)
# [1] 0.33 0.42 0.33 0.38 0.47 0.50 0.71 0.46 0.41 0.55 0.49 0.31 0.30 0.43 0.25

as.numeric(Achange)            # [1] 1 0 1 1 1 1 0 0 1 0 0 0 1 0     (8 runs)
runs.test(as.factor(Achange))

Runs Test

data:  as.factor(Achange)
Standard Normal = 0, p-value = 1
alternative hypothesis: two.sided

as.numeric(Ahigh)             # [1] 0 0 0 0 1 1 1 1 0 1 1 0 0 1 0     (7 runs)
runs.test(as.factor(Ahigh))

Runs Test

data:  as.factor(Ahigh)
Standard Normal = -0.7898, p-value = 0.4297
alternative hypothesis: two.sided


So group A looks random according to both methods. (Bear in mind that the test relies on large sample theory, and your N is 15, but the data seem OK. Also remember that, if you have a bunch of groups and you run a bunch of tests on them, something is likely to show up whether its real or not.)

2) To determine if the distributions are similar, you can calculate summary statistics and plot graphs. I often find graphs most helpful. Two graphs that I like are kernal density plots and qq-plots. A kernal desity plot is like a smoothed histogram. It's easy to plot several distributions together. A qq-plot is a scatterplot of one distribution against another after both have been sorted into ascending order. If the two distributions have similar shapes, the points should lie along a 45 degree line through the origin. Most people think of qq-plots as a way to compare a distribution to a theoretical normal distribution, but they can be used for any theoretical distribution (e.g., Weibull) or the empirical distribution of another data set. Thus, you could make a series of plots for pairwise comparisons. It also helps to plot the 45 degree line to help with interpretation. Here is some R code and graphs:

summary(Aprop)
#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
# 0.2456  0.3281  0.4155  0.4215  0.4805  0.7111

windows()
plot( density(Aprop), col="gold", xlim=c(0,1))
lines(density(Bprop), col="red")
lines(density(Cprop), col="blue")
legend("topleft", c("Aprop","Bprop","Cprop"), lty=1,
col=c("gold","red","blue"))

windows()
qqplot(Aprop, Cprop)
abline(0, 1)

windows()
qqplot(Aprop, Bprop)
abline(0, 1)


The first qq-plot shows that the distributions of the A group and the C group are pretty similar--they follow the 45 degree line pretty well. It's just that C is consistently higher than A. This can be confirmed by looking at the density plot. At first, the lower qq-plot looks OK, but you notice that the 45 degree line doesn't even show up in the plot window. That's a big tip-off. If you look at the scales, you see that B varies between .91 and 1.0, whereas A varies from .3 to .7. A and B are just not that similar, as you can see in the density plot. Density plots are really informative, but they can also exaggerate the shape of a distribution. That's why it's nice to plot both. A and C look rather different in the density plot, but based on the qq-plot, I'd bet the differences in the shapes of the distributions wouldn't be reliable. And both plots suggest @Zach is right about that outlier.

-

i'm not sure if i understand the problem correctly - so correct me if i'm wrong but here's my thoughts:

So, is spottiness non-random for members in a group? And does spottiness vary significantly between the groups?

I'm not sure if the first question is answerable with the data presented. If the data was available in the following form:

         Head Tail Body Group
Lizard 1   Y    N    N    A
Lizard 2   Y    Y    N    A
...


One would be able to test the randomness of individuals. However, the individuals are - as far as I understand not available - right? The same applies for the second question which I would simply test with a Mann-Whitney-U test or something like that (in R use ?wilcox.test).

So I hope i got your problem right - if not please correct me!

EDIT

As for the second question: whether to determine if you have significantly more individuals with spots in one group than in another. one (maybe to simple) way was taking the distributions given above (with normal densities in blue):

and assuming a normal distribution and then testing each group against the mean. that would be: $z = (A_1 - \bar{A}) / \sigma_A$ with $\bar{A}$ being the mean of group $A$.

Correct me if I made a mistake but this would give the following p-values:

pnorm((A-mean(A))/sd(A))
[1] 0.21718755 0.47952112 0.20871568 0.35415346 0.66053450 0.74852194
[7] 0.99325922 0.61952780 0.47553465 0.85819411 0.72317286 0.16977245
[13] 0.14707919 0.52412299 0.06677287

-
Hi - the data is also available in the form Lizard 1,Y,N,N,GroupA, Lizard2,Y,Y,N,GroupB <-- Ie presence/absence only for each lizard. Correct me if I'm wrong, but if presence/absence were random, then if GroupX had 100 lizards, we would predict a distribution of 50,50,50... So is there not some way to see whether an observed distribution of 30,60,60 is significant? –  Trees4theForest Dec 18 '11 at 3:15
if you had the vector for each group in the above form - you may test wilcox.test(A,B) in order to see wether the two groups have different location parameters... was this your question? –  Seb Dec 18 '11 at 11:58
i edited my answer a bit - hope this helps. –  Seb Dec 18 '11 at 12:38