I’m somewhat lost as to how to approach a problem I have and was hoping someone can suggest the most appropriate method.
I have a data set consisting of concentration values for 8 molecules. These are measured on approx. 230 samples from approx. 70 patients. For many patients there are multiple samples taken at different dates. However, this isn’t a controlled timecourse but independent (presumably, but probably not quite) events for each patient. As such the number and time of samples for each patient differ significantly so this data would not fit into a nice balanced repeated measures model.
For each sample we can assign 2 clinical phenotypes, let’s call them A and B, which are presumably independent (although probably not entirely). So we could assign 4 categories based on A or B, but we could also sub classify both A and B based on severity. For now, I think binary categories would be sufficient. Some preliminary analysis someone else did on the data suggests that these molecules will not be able to distinguish A, but may be able to distinguish B, independent of A. An additional complication is that patients come from 2 different institutes (about half from each). Also, the data appear to deviate fairly significantly from normality, even after a log transformation. Most of the variances (comparing B+ vs B-) are similar after log transformation.
The goal is to determine if any of the 8 molecules, or some combination of them, can predict clinical phenotype B. Ideally, predict A and B, but B is the more important one for this study.
My initial thought is to do some kind of mixed model analysis, two-way across A and B, with institute and patient as random effects. This would help determine if A and B are truly independent and if there is any significant confounding effect from institute. This should probably use a non-parametric test.
So does this make sense (my confusion suggests otherwise), and if so, what exactly is the test I should be using (I’m using R for most of my analysis)? If not, what is a better way to look at this data?
Response to Chi's questions:
When I plot the molecule concentrations they are quite skewed (even after a log2 transformation) and Shapiro-Wilk Normality Tests on each molecule (B- vs B+ samples) fails to support normal distribution. Therefore, I felt parametric tests would not be justified due to deviation from an assumed normal distribution.
On the other hand, I've read that with sample sizes > 30, due to central limit theory (something I'll have to refresh myself on...) parametric tests are still applicable even if the underlying distribution is non-normal. In my case there are around 170 B- and 60 B+ samples in total (although they are not from individual patients) so perhaps parametric tests would be the best choice? This is something I'm still confused about.
I did do, as a preliminary test, both t-test and Mann-Whitney tests between B- and B+ for each molecule, followed by a Benjamini & Hochberg FDR adjustment (Which raises another question: is B & H appropriate for such a small number of tests?). With a cutoff of q < 0.05, I got a significant difference for 3 molecules from the Mann-Whitney test, whereas only 2 of those were sig from the t-test of log2 transformed data. Not certain what this means. Thanks!
I have noticed there are also some identical concentration across several samples for some molecules. I don't fully understand the source of this artifact, but was told had something to do with data validation. You might see this with readings that are below the sensitivity of the assay, but I'm seeing identical values that are not all the lowest value. On a Q-Q plot these show up as horizontal lines and I'm not sure how this will affect the analysis.
I tried to set up a linear mixed effects model in R but just could not get it to work. Seems like the 2 most used packages are nlme or lme4. I can't figure out the syntax of the formulas and so far despite lots of searching haven't found a thorough explanation that I understand.
Lets call my variables phen.a and phen.b for the independent variables A and B, above, instit and patient for institute and patient id. instit and patient are factors, phen.a dn phen.b are logical but could be changed to factors, especially we decide to add severity grades.
What I think I should be doing is a mixed effects model looking at: fixed effects: phen.a, phen.b, interaction:phen.a x phen.b random effects: instit, patient
If someone could provide the formulas I need that would be extremely helpful.
In the meantime, what I've tried, since I could get it to work, is an ANOVA looking at the 2 independent variables, their interactions, and additional interactions to phen.b, which is of particular interest. The formula I used in the linear model was:
concentration ~ phen.a*phen.b + phen.b:instit + phen.b:patient
where concentration is the log2 of the concentration -- repeated for each of the 8 molecules. I then did a p.adjust(..., method="FDR") on each of the 5 resulting p values to correct for the 8 comparisons.
This resulted in 1 molecule being quite significant for phen.b (with no sig for phen.a and no interactions), and a second molecule sig for phen.a, with an interaction between phen.b:instit. Not sure how to interpret this interaction on a non-significant independent variable.
So now another question is whether this model is valid and if the interactions I chose make sense. I'm also wondering since we are essentially looking at 2 x 2 categories if linear models are appropriate or if there is some kind of chi2 type test that would be better (although the dependent variable is continuous).
Another update: So I managed to learn how to use the lme function and tried the model:
fixed = concentration ~ phen.a * phen.b * instit
random = ~1 | patient
to find any effect from phen.a or phen.b and determine if such effects are confounded by instit, using the random patient effect to account for differences in intra vs inter patient variation.
However, as I played with the data I noticed that while the concentration is a continuous variable there are 1 or 2 particular values that show up for many of the samples, sometimes accounting for half or more of all samples, and completely skewing the distribution. It turns out that these values actually represent the minimal detectable concentrations by the assays (in some cases 2 assays were used, accounting for the 2 different values that show up repeatedly). The way the data were coded, the minimal detectable concentration is used for any sample for which that molecule could not be detected. So that value really represents: <= x.
How do I deal with this type of minimal cutoff? Is it fair to leave them as is, for which I found a couple of significant results?
I tried just excluding these values but then I lose all significant results. I don't think that's valid, though, because the values aren't missing, they really represent a molecule that is either not present or present at very small concentrations. If there is a predominance of very low concentration correlated with a phenotype that is a very important observation.