If you want to know whether treatments and phenotypes are independent
categorical variables, then consider the counts (exclusive of totals)
to make a contingency table, and do a chi-squared test of independence.
First, let's temporarily ignore the double counting and pretend that the
number of subjects per row is the total of the counts per row. That amounts
to pretending there were 840 subjects.
TABLE = matrix(c(396, 44, 10, 54, 67, 37, 0, 13,
62, 19, 8, 22, 45, 13, 0, 5,
40, 5, 0, 0), byrow=T, nrow = 5)
TABLE
[,1] [,2] [,3] [,4]
[1,] 396 44 10 54
[2,] 67 37 0 13
[3,] 62 19 8 22
[4,] 45 13 0 5
[5,] 40 5 0 0
sum(TABLE)
[1] 840
The expected counts are $E_{ij} = R_iC_j/G,$ where $R_i, C_j, G$ are
corresponding row and column totals, and the grand total $G=840.$ Then under $H_0$ that row and column categories are independent, the chi-squared statistic is
$$Q = \sum_{i,j}\frac{(X_{ij}-E_{ij})^2}{E_{ij}} \stackrel{\text{aprx}}{\sim}
\mathsf{Chisq}(\nu = (r-1)(c-1)),$$
where the approximation depends of having "most" $E_{ij} > 5$ and all
of them larger than 3. Here $\nu = (5-1)(4-1) = 12.$
However, for our data, several of the $E_{ij}$ in the lower-right corner
are too small. In this case upon request, R will simulate the distribution
of $Q$ to give a reasonably accurate P-value. Omitting the false start with
a warning about $E_{ij}$'s being to small and moving directly to simulation,
we have a highly significant result with P-value near 0 (as small as the
number of iterations of the simulation permits). Without the simulation
the P-value is also very near 0, 1.478e-12.
[Notice that $\nu$ is missing from the output because the chi-squared distribution was not used.]
Results = chisq.test(TABLE, simulate.p=T); Results
Pearson's Chi-squared test with simulated p-value (based on 1e+05 replicates)
data: TABLE
X-squared = 82.334, df = NA, p-value = 1e-05
However, this result is questionable because of the multiple counting, so we
need to try to adjust for multiple counting in some way.
Below I propose a few possible methods of adjustment. I have not encountered
this problem before, so maybe someone can find a better method than any of them.
Adjust statistic Q: Suppose for the moment that every count $X_{i,j}$
were doubled. Then every $E_{i,j}$ is doubled and so is the statistic $Q.$
Thus, a 'cure' for this systematic doubling of counts would be to halve $Q$
before finding the P-value based on $\nu$ degrees of freedom. In our case according to your Comment,
we have sporadic multiple counts that we might view as inflating $Q$ by
a factor of $c = 840/701 \approx 1.2.$ We could use an adjusted statistic
$Q^\prime = Q/c = 82.334/1.12 = 73.51,$ and the corresponding P-value would again be very nearly 0 (7.02 e-11).
1 - pchisq(73.51, 12)
[1] 7.021073e-11
An objection to this method of adjustment is that our multiple counts are not
uniformly spread across the table. For example, there are no multiple counts
for Treatment 5, while Treatment 1 has $504 - 434 = 70$ multiple counts.
Its strong point is that it does not interfere with overall proportions of counts.
Adjust treatment counts by "exclusion." Sometimes when crucial data values are missing,
one uses 'imputation' to estimate them. Here we have the opposite problem, we
have unwanted extra counts. We might seek a rational way to exclude them.
rowSums(TABLE)
[1] 504 117 111 63 45
Scaling method: A simple way to make an adjustment would be just to scale each row. However, due to rounding, the adjusted row total may be slightly incorrect
(below for Treatment 1, 436 instead of 435).
t1.over = TABLE[1,]; t1.over; sum(t1.over)
[1] 396 44 10 54
[1] 504
t1.adj = round(t1.over*435/504); t1.adj; sum(t1.adj)
[1] 342 38 9 47
[1] 436
Then re-scale other the other rows with overcounts.
Random method: Alternatively, we might imagine original researchers looking at an individual in Treatment 1 and
trying to decide whether to classify it as phenotype A or B. Unable to decide
for sure, they just count it as both. So we are left with counts
t1.over = c(369, 44, 10, 54) totaling 504, while the number of subjects
in Treatment 1 is known to be 435.
We could use the sample function in R
to 'decide' for them on a probability basis.
t1.s = sample(1:4, 435, repl=T, prob=TABLE[1,])
table(t1.s)
t1.s
1 2 3 4
346 40 8 41
t1.radj = as.numeric(as.matrix(table(t1.s))); t1.radj
[1] 346 40 8 41
sum(t1.radj)
[1] 435
This random method of exclusion might be done repeatedly in a simulation to see
how much P-values vary.
