I have the allele frequencies of both case and control groups, and I want to see if the allele frequencies are significantly different or not, I want to do the analysis in R.
You simply have 2 groups and want to compare frequencies?
(1) BY NORMAL APPROXIMATION
For n*p > 10 and n*(1-p) > 10 the normal approximation works well. The you can calculate the mean difference in the estimated proportions and devide this by its standard error
SE(p1-p2) = sqrt(p1*(1-p1)/n1 + p2*(1-p2)/n2)
to get the z-value. The confidence interval and the test can be done using the normal distribution (function qnorm(.) for the CI and pnorm(.) to get a p-value for the test).
(2) BY FISHERS EXACT TEST
The exact way to compare the proportions can be achieved using Fishers exact test:
fisher.test(.)
The odds ratio from Fisher's exact test can be transformed to a relative risk by
RR = OR / (1 – p + (p*OR))
where OR is the odds ratio, RR the relative risk and p the proportion in the control group.
(3) BY A LOGISTIC REGRESSION
If you have a more complicated design with some covariables, a logistic regression (ie. a generalized linear model of the binomial family with logit link) would be most appropriate. Then R function is glm(.)
You simply have 2 groups and want to compare frequencies?
(1) BY NORMAL APPROXIMATION
For n*p > 10 and n*(1-p) > 10 the normal approximation works well. The you can calculate the mean difference in the estimated proportions and devide this by its standard error
SE(p1-p2) = sqrt(p1*(1-p1)/n1 + p2*(1-p2)/n2)
to get the z-value. The confidence interval and the test can be done using the normal distribution (function qnorm(.) for the CI and pnorm(.) to get a p-value for the test).
(2) BY FISHERS EXACT TEST
The exact way to compare the proportions can be achieved using Fishers exact test:
fisher.test(.)
The odds ratio from Fisher's exact test can be transformed to a relative risk by
RR = OR / (1 – p + (p*OR))
where OR is the odds ratio, RR the relative risk and p the proportion in the control group.
(3) BY A LOGISTIC REGRESSION
If you have a more complicated design with some covariables, a logistic regression (ie. a generalized linear model of the binomial family with logit link) would be most appropriate. Then R function is glm(.)
Thank you Jochen Wilhelm , my data is in an excel data sheet consist of 2 columns, one the details of genotype from case group and the other control group.( The data is like this, A/A, A/G, A/A, G/G ...........)
Now how can analyze the data in R to detect the genotype distribution difference (if any) between these 2 groups.
Given the columns in Excel have the headers "cases", and "controls", select the range including the headers and copy it it the clipboard.
Read the data in R with
x = read.delim("clipboard")
transform the data into a more convenient format (one column per variable, one row per observation):
x = reshape(x,dir="long",varying=1:2,v.names="genotype",times=colnames(x),timevar="group")
get the contingency table with
tab = xtabs(~group+genotype,x)
this shows you the distribution. You may visualize it with a mosaic plot:
mosaicplot(tab)
You can test the hypothesis if both groups are samples from the same population by
fisher.test(tab)
Note that you dont get an odds ratio when you don't have a 2x2 table.
To perform a logistic regression (e.g. to get estimates for log odds ratios for the genotypes) you may google for
STATS 330: Lecture 28 Two-Dimensional Contingency Tables
the main problem will be to set up a data.frame in a form that contains variables and counts as required by the glm(.)-function:
df = as.data.frame.table(tab)
so you can get the model by
model = glm(Freq~group*genotype,df,family=poisson)
The rest and the interpretation of the coefficients and all that can be found in the "STATS 330" presentation.
(Note that here we dela with a Poisson-model as we have the allele *frequencies* (= counts = Poisson distribution) from the contingency table)
Thank you very much Jochen Wilhelm, for your explanation, I really appreciate it. I put all genotypes data from case and control groups in a single column (header=all) in excel, and in another column (header=dis), I entered "0" for each genotype from control group and "1" for each genotype from case group. Then I imported it to R and analyzed with "SNPassoc" package wirh "association(dis~all, data=mydat)" command. Did I do the right thing?
- Badges
- Science topic
- Similar topics
- Mathematics
- Statistics
- Statistical Software
- R Statistical Package