I am conducting a mulitple first order regression analysis of genetic data. The vectors of y-values do not all follow a normal distribution, therefore I need to implement a non-parametric regression using ranks.

Is the lm() function in R suitable for this, i.e.,

lin.reg <- lm(Y~X*Z)

where Y, X and Z are vectors of ordinal categorical variables?

I am interested in the p-value assigned to the coefficient of the interaction term in the first order model. The lm() function obtains this from a t-test, i.e., is the interaction coefficient significantly different from zero.

Is the automatic implementation of a t-test to determine this p-value appropriate when the regression model is carried out on data as described?



Sample data for clarity:

Y <- c(4, 1, 2, 3) # A vector of ranks
X <- c(0, 2, 1, 1) # A vector of genotypes (0 = aa, 1 = ab, 2 = bb)
Z <- c(2, 2, 1, 0)
  • $\begingroup$ Could you explain what are $Y$, $X$, and $Z$? $\endgroup$ – chl Jan 19 '11 at 20:51
  • $\begingroup$ Y is a vector of gene expression values. These are not gaussian, thus they are represented as a vector of ranks. X and Z are vectors of SNP data (genotypes at a particular location) for multiple samples. They are encoded as 0, 1, and 2 representing aa, ab and bb respectively. $\endgroup$ – Darren J. Fitzpatrick Jan 19 '11 at 20:54
  • $\begingroup$ Is that 'first order regression' in the sense of Karalič and Bratko (1997)?? dx.doi.org/10.1023/A:1007365207130 If not, could you explain what you do mean by 'multiple first order regression analysis'? $\endgroup$ – onestop Jan 19 '11 at 20:55
  • $\begingroup$ No, it is first order in the sense that Y ~ X + Z + XZ + e, where X and Z are main effects, AB is the effect of the interaction and e is the error. Beta 1, 2 and 3 coefficients are omitted for clarity. $\endgroup$ – Darren J. Fitzpatrick Jan 19 '11 at 20:59
  • 1
    $\begingroup$ 1) "The vectors of y-values do not all follow a normal distribution, therefore I need to implement a non-parametric regression" is a non-sequitur. Linear regression is reasonably rubust to departures from normality with large samples, and as you're in genetics I suspect you have a very large sample. And non-parametric methods are not the only alternative; you could apply a transformation, or use generalized linear models. 2) You still haven't given us any idea information about $Z$, e.g. how many levels does it have? $\endgroup$ – onestop Jan 19 '11 at 21:56

If your response variable is ordinal, you may want to consider and "ordered logistic regression". This is basically where you model the cumulative probabilities {in the simple example, you would model $Pr(Y\leq 1),Pr(Y\leq 2),Pr(Y\leq 3)$}. This incorporates the ordering of the response into the model, without the need for an arbitrary assumption which transforms the ordered response into a numerical one (although having said that, this can be a useful first step in exploratory analysis, or in selecting which $X$ and $Z$ variables are not necessary)

There is a way that you can get the glm() function in R to give you the MLE's for this model (other wise you would need to write your own algorithm to get the MLEs). You define a new set of variables, say $W$, where these are defined as

$$W_{1jk} = \frac{Y_{1jk}}{\sum_{i=1}^{i=I} Y_{ijk}}$$ $$W_{2jk} = \frac{Y_{2jk}}{\sum_{i=2}^{i=I} Y_{ijk}}$$ $$...$$ $$W_{I-1,jk} = \frac{Y_{I-1,jk}}{\sum_{i=I-1}^{i=R} Y_{ijk}}$$

Where $i=1,..,I$ indexes the $Y$ categories, $j=1,..,J$ indexes the $X$ categories, and $k=1,..,K$ indexes the $Z$ categories. Then fit a glm() of W on X and Z using the complimentary log-log link function. Denoting $\theta_{ijk}=Pr(Y_{ijk}\leq i)$ as the cumulative probability, the MLE's of the theta's (assuming a multi-nomial distribution for $Y_{ijk}$ values) is then

$$\hat{\theta}_{ijk}=\hat{W}_{ijk}+\hat{\theta}_{(i-1)jk}(1-\hat{W}_{ijk}) \ \ \ i=1,\dots ,I-1$$

Where $\hat{\theta}_{0jk}=0$ and $\hat{\theta}_{Ijk}=1$ and $\hat{W}_{ijk}$ are the fitted values from the glm.

You can then use the deviance table (use the anova() function on the glm object) to assess the significance of the regressor variables.

EDIT: one thing I forgot to mention in my original answer was that in the glm() function, you need to specify weights when fitting the model to $W$, which are equal to the denominators in the respective fractions defining each $W$.

You could also try a Bayesian approach, but you would most likely need to use sampling techniques to get your posterior, and using the multinomial likelihood (but parameterised with respect to $\theta_{ijk}$, so the likelihood function will have differences of the form $\theta_{ijk}-\theta_{i-1,jk}$), the MLE's are a good "first crack" at genuinely fitting the model, and give an approximate Bayesian solution (as you may have noticed, I prefer Bayesian inference)

This method is in my lecture notes, so I'm not really sure how to reference it (there are no references given in the notes) apart from what I've just said.

Just another note, I won't harp on it, but I p-values are not all they are cracked up to be. A good post discussing this can be found here. I like Harlod Jeffrey's quote above p-values (from his book probability theory) "A null hypothesis may be rejected because it did not predict something that was not observed" (this is because p-values ask for the probability of events more extreme than what was observed).

  • $\begingroup$ @probabilityislogic It's not clear to me whether your model treats all variables as categorical ones or not, and what is the advantage of modeling ranks using an adjacent or cumulative-logit model. Anyway, when dealing with SNP data, we often assume an additive model (aka, allelic dosage) and treat the frequency of minor allele {0,1,2} as a continuous predictor. $\endgroup$ – chl Jan 20 '11 at 9:13
  • $\begingroup$ The questioner's sample data has as many values for Y as there are observations. If this is true of his real data, which I guess it is as they are ranks (ignoring any ties), maximum likelihood estimation isn't going to work, as the number of $\theta$s to estimate increases with the numbers of observations. $\endgroup$ – onestop Jan 20 '11 at 13:28
  • $\begingroup$ @chl, This model treats X and Z as categorical in the standard ANOVA-regression (factor contrasts) representation for categorical regressors. That's why I have the $j$ and $k$ subscripts in the equations (if they were continuous only the $i$ subscript would be there). And from what it sounds like your "minor allele" you have information about the scale of the categories, the "ordinality" of X and Z only affect the interpretation of coefficients. But by treating it as a numerical value, you are effectively replacing "degrees of freedom" with "model assumptions". $\endgroup$ – probabilityislogic Jan 21 '11 at 16:56
  • $\begingroup$ @onestop - if this is true ($n=I$), then the model as I have described it is unidentifiable which means that the problem basically needs more structure in order to get a unique solution. One way to add this "structure" is through the Bayesian approach, setting a slightly informative prior on the probabilities (such as the uniform prior). $\endgroup$ – probabilityislogic Jan 21 '11 at 17:07
  • $\begingroup$ @probabilityislogic Thanks for clarifications. If you're interested in how to "scale" SNP data, there was a nice paper by Waaijenborg & Zwinderman with a pretty nice use of optimal scaling, Correlating multiple SNPs and multiple disease phenotypes: penalized non-linear canonical correlation analysis (Bioinformatics 2009 25(21):2764). I'm still not clear myself with how we could convert a rank outcome to an ordinal one (hence echoing @onestop's comment), so I'm waiting for an helpful OP's feedback. $\endgroup$ – chl Jan 22 '11 at 11:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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