1
$\begingroup$

I'm working on a project with genetics but I think my problem is applicable to general statistics.

I want to test frequency (Minor Allele Frequency) of a SNP/variant across 5 age categories to see if there is any change. I am expecting either linear or a 'U' shaped curve if such an association exists. I have about 1000 individuals within each age group and have calculated the allele frequency for each age group (ranging from 0.001-0.5), hence I only have 5 frequencies in total.

Now the problem is that I have 2 million variants/SNPs and so I can't plot each SNP to see the trend before I commence the regression.

Currently I propose a linear model with Frequency ~ Age to see if there is a linear trend and Frequency ~ Age + I(Age^2) to see if there is a quadratic curve. Then using a likelihood ratio test to see whether the quadratic term is significant and deducing whether there is a curve.

My question is whether this is the best approach, statistically speaking?

$\endgroup$
8
  • $\begingroup$ Erika, what does a 5x10^-8 threshold mean? You want the p-value to be that value? That would mean that your relationships needs to be nearly exactly like a polynomial relation. For instance the example from @olooney with 250 200 200 220 250 has "only" a p-value of 0.001. You will be shooting your shotgun with 2 million pieces hoping to find a near-perfect linear or polynomial curve, but that "luck" may not be because you found an actual hit. Your test has so little power that you need extremely many trials to be lucky.... $\endgroup$ Dec 18, 2018 at 9:44
  • $\begingroup$ ... Also the interpretation of the p-value is strange. What if you find a weak linear relationship like 250 251 252 253 254 but with the perfect p-value? What if you find negative coefficients for the quadratic relationship (e.g. inverted U or J shapes). What is a U shape anyway? Why must it be a quadratic relationship? Any deviation from a quadratic relationship will give you p-values that are not passing the threshold standard of your field. $\endgroup$ Dec 18, 2018 at 9:52
  • $\begingroup$ @MartijnWeterings The 5x10^-8 threshold is based on the Bonferroni correction where only p-values less than this are deemed significant, to control the number of false positives when performing so many tests. We hypothesize that certain SNPs have a 'U' shape in frequency across age, i.e. monotonically decease with age and then increase, and certain SNPs that increase/decrease with age. Our goal is to identify these SNPs. Do you have any suggestions of a better analysis method? I've considered increasing the number of age bins in the hope of increasing power. $\endgroup$ Dec 18, 2018 at 17:14
  • $\begingroup$ There are possibly many better methods since looking for a p-value needle in the haystack that is modeled by a polynomial (thus excluding needles from different models and needles not found) may give you just a straw and mis many needles. But, for better methods, you need to clearly defined what you are looking for and why you are looking for this. That you arw just looking for a linear relationship (which can be a flat line as well) or a polynomial relationship is hard to believe. $\endgroup$ Dec 18, 2018 at 17:22
  • $\begingroup$ @MartijnWeterings I understand. I am investigating the relationship between frequency & age, for 2 million SNPs. I am ultimately looking to identify SNPs where frequency decreases with age until middle age and then begins to increase, hence the 'U' shape on a plot of Frequency~Age. I am also trying to identify SNPs where frequency either increases or decreases with age, hence the linear relationship. My understanding was that a flat line would not produce a significant p-value for the linear model vs. null as this would show age does not have a significant effect on Frequency $\endgroup$ Dec 18, 2018 at 21:33

1 Answer 1

2
+50
$\begingroup$

You can do somewhat better by using logistic regression instead of linear regression, but otherwise your approach is fine.

Note that with the linear model you propose, the model has no way of "knowing" that there are 1000 observations in each age band because it only "sees" the average. That means the linear model is forced to estimate variance purely from the 5 provided data points which is not very reliable. (It's also extremely likely that we will be violated the underlying assumption of normally distributed error term, although it would of course be extremely difficult to tell one way or the other with just 5 data points!) In contrast to these difficulties, a logistic regression model will be able to "see" the sample size and exploit the known relationship between n, p, and variance of the binomial distribution and we won't have to worry about violating the normality assumptions of the linear model.

To adapt your approach to logistic regression, you'll need to use a categorical variable (a factor) as the response variable instead of a continuous frequency. One way to do this is simply to replicate each row of data the appropriate number of times but another faster way to do this is with observation weights.

Let's say we format the data to look like this:

| age | variant | n_people | 
|-----|---------|----------| 
| 20  | 0       | 750      | 
| 20  | 1       | 250      | 
| 30  | 0       | 800      | 
| 30  | 1       | 200      | 
| 40  | 0       | 800      | 
| 40  | 1       | 200      | 
| 50  | 0       | 780      | 
| 50  | 1       | 220      | 
| 60  | 0       | 750      | 
| 60  | 1       | 250      | 

Then (using R, for example) we can build models which include the quadratic age term, the main (linear) effect only, and the reduced model with only an intercept.

full_model <- glm(variant ~ age + I(age^2), data=df, weights=n_people, family=binomial())
main_effects_model <- glm(variant ~ age , data=df, weights=n_people, family=binomial())
reduced_model <- glm(variant ~ 1, data=df, weights=n_people, family=binomial())

The full model looks like so (using the fake data I mocked up):

> summary(full_model)
Call:
glm(formula = variant ~ age + I(age^2), family = binomial(), 
    data = df, weights = n_people)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-20.960  -19.353    3.096   25.931   26.491  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.2617126  0.2978734  -0.879 0.379616    
age         -0.0576910  0.0163196  -3.535 0.000408 ***
I(age^2)     0.0007348  0.0002017   3.643 0.000270 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 5319.2  on 9  degrees of freedom
Residual deviance: 5305.8  on 7  degrees of freedom
AIC: 5311.8

Number of Fisher Scoring iterations: 5

Note that R is helpfully performing a z-test for each coefficient, which already suggests that the quadratic age term may be statistically significant (again, for the fake data.) However, these z-tests on individual coefficients are an approximation at best, so we should use the likelihood ratio test, which is the gold standard.

Here's an example in R of using the LR test to compare the full model to the reduced model. This tells us if including age (both linear and quadratic terms) is better than excluding age from the model altogether:

> anova(full_model, reduced_model, test="LR")
Analysis of Deviance Table

Model 1: variant ~ age + I(age^2)
Model 2: variant ~ 1
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)   
1         7     5305.8                        
2         9     5319.2 -2  -13.472 0.001188 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Alternatively, we could compare the full model to the main effects model with anova(full_model, main_effects_model) if we are specifically interested whether or not the inclusion of the quadratic term to the model is statistically significant:

> anova(full_model, main_effects_model, test="LR")
Analysis of Deviance Table

Model 1: variant ~ age + I(age^2)
Model 2: variant ~ age
  Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
1         7     5305.8                          
2         8     5319.0 -1  -13.242 0.0002738 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Note that there are many approaches to this same problem, say a 5x2 $\chi^2$ test, but the logistic regression approach I've outlined above seems closest to your original proposal while avoiding some of the theoretical concerns of the linear model.

$\endgroup$
2
  • $\begingroup$ This is perfect, thank you for explaining so clearly. As I mentioned I will be performing this regression for 2 million different SNPs and so will need to apply a multiple testing correction to the p-values I obtain (the standard in this field is a 5x10^-8 threshold). The objective is to identify SNPs that have significant linear or U shape trends. Am I therefore correct in thinking the p-value I should take is that obtained from the likelihood ratio test rather than the coefficients in the model? $\endgroup$ Dec 18, 2018 at 0:17
  • $\begingroup$ Definitely. This answer stats.stackexchange.com/a/48683/48250 goes into more detail. $\endgroup$
    – olooney
    Dec 18, 2018 at 14:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.