-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathProjectCode.R
490 lines (370 loc) · 21.3 KB
/
ProjectCode.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
# Merge the genotype & neuro data ----
#getwd()
#setwd("C:/Users/RA/Desktop/Aanya Data") <- for windows!
#setwd for mac by going to Session/Set Working Directory/Choose Directory
#On a mac, press control + L to clear the console
getwd()
#setwd("C:/Users/RA/Desktop/Aanya Data")
setwd("/Volumes/EC_RAs/Aanya - COMTBDNF/Aanya Data/")
dir()
genotype_data <- read.csv("Genotype Data.csv")
neuro_data <- read.csv("SOLVED Neuro Data January 22 2020.csv")
dataset <- merge(genotype_data, neuro_data, by="ID")
# Install packages -----
install.packages("carData")
install.packages ("car")
install.packages("dplyr")
install.packages("tidyverse")
install.packages("FSA")
install.packages("ggpubr")
install.packages("HardyWeinberg")
install.packages("afex")
install.packages("knitr")
install.packages("arsenal")
install.packages("magrittr")
# Load packages ----
library(carData)
library(car)
library(dplyr)
library(tidyverse)
library(FSA)
library(emmeans)
library(ggpubr)
library(HardyWeinberg)
library(afex)
library(knitr)
library(arsenal)
library(magrittr)
# Assigning Groups ----
dataset$group_membership[dataset$Group == 1 & dataset$HRTever == 0] <- "BSO"
dataset$group_membership[dataset$Group == 1 & dataset$HRTever == 1 & dataset$E2now == 1] <- "BSO-E2"
dataset$group_membership[dataset$Group == 3 & dataset$HRTever == 0] <- "AMC"
dataset$group_membership[dataset$Group == 1 & dataset$HRTever == 1 & dataset$HRTnow == 0] <- "BSO_HRT_Past"
dataset_new <- subset(dataset, !group_membership=="BSO_HRT_Past")
dataset_new <- dataset_new[!is.na(dataset_new$group_membership),]
dataset_new <- dataset_new[!is.na(dataset_new$ID),]
dataset_new <- dataset_new[!is.na(dataset_new$SPWM_WME_T1),]
#View(dataset_new)
#Set white spaces into NAs
#dataset_new[dataset_new == ""] <- NA
#View(dataset_new)
# Adding BDNF Met Carriers to the dataset ----
dataset_new$BDNF_combined[dataset_new$BDNF == "Val"] <- "Val"
dataset_new$BDNF_combined[dataset_new$BDNF == "Het" | dataset_new$BDNF == "Met"] <- "Met_carrier"
#use clean_data as a dataset
clean_data <- read.csv("clean_data.csv")
clean_data <- dataset_new
# Group Counts for demographics table ----
# Use the %>% (pipe operator) to simplify code which can otherwise have a lot of parentheses!
clean_data %>% group_by(group_membership) %>% tally()
clean_data %>% group_by(group_membership, COMT) %>% tally()
# get a tally for BDNF numbers
clean_data %>% group_by(group_membership, BDNF_combined) %>% tally()
# get a tally for the APOE #s
clean_data %>% group_by(group_membership, e4) %>% tally()
#get a tally for the PdG numbers
clean_data %>% group_by(group_membership, PdG) %>% tally()
#get a tally for the cancner numbers
clean_data %>% group_by(group_membership, Cancer_Treatments) %>% tally()
# Standard error = SD / sqrt(n-1)
# Introducing the sem calculation into the summarize funtion makes your life a lot easier - eliminates many lines of code!
View(clean_data)
clean_data %>% group_by(group_membership) %>% summarize(mean=mean(E1G, na.rm = T),
sd=sd(E1G, na.rm = T),
n=length(Group),
sem=sd(E1G, na.rm = T)/sqrt(n-1))
clean_data %>% group_by(group_membership) %>% summarize(mean=mean(MenopauseAge, na.rm = T),
sd=sd(MenopauseAge, na.rm = T),
n=length(Group),
sem=sd(MenopauseAge, na.rm = T)/sqrt(n-1))
clean_data %>% group_by(group_membership) %>% summarize(mean = mean(EduYears, na.rm = T),
sd =sd(EduYears, na.rm = T),
n=length(Group),
sem=sd(EduYears, na.rm = T)/sqrt(n-1))
clean_data %>% group_by(group_membership) %>% summarize(mean=mean(Age, na.rm = T),
sd=sd(Age, na.rm = T),
n=length(Group),
sem=sd(Age, na.rm = T)/sqrt(n-1))
clean_data %>% group_by(group_membership) %>% summarize(mean=mean(BMI, na.rm = T),
sd=sd(BMI, na.rm = T),
n=length(Group),
sem=sd(BMI, na.rm = T)/sqrt(n-1))
clean_data %>% group_by(group_membership) %>% summarize(mean(Cancer_Treatments, na.rm = T), sd(Cancer_Treatments, na.rm = T))
clean_data %>% group_by(group_membership) %>% summarize(mean=mean(PdG, na.rm = T),
sd=sd(PdG, na.rm = T),
n=length(Group),
sem=sd(PdG, na.rm = T)/sqrt(n-1))
# Group comparison tests for demo ----
# One factorial ANOVA
#Age
attach(clean_data)
Age.aov<-aov(Age~group_membership)
summary(Age.aov)
leveneTest(Age~group_membership) #assess the equality of the variances
# Extract the residuals
Age_aovresiduals <- residuals(object = Age.aov )
shapiro.test(x = Age_aovresiduals )
#Does pass Levene's Test and Shapiro
#No. With History of Cancer Treatment
Cancer_Treatment.aov <- aov(Cancer_Treatments~group_membership, data = clean_data)
summary(Cancer_Treatment.aov)
leveneTest(Cancer_Treatments~group_membership) #Assess the quality of the variances
#Extract the residuals
Cancer_Treatments_aovresiduals <- residuals(object = Cancer_Treatment.aov)
shapiro.test(x = Cancer_Treatments_aovresiduals)
#Is significant so run a Kurskal-Wallis Rank Sum Test
kruskal.test(Cancer_Treatments~group_membership)
# Post-hoc for KW is the Dunn Test
clean_data$group_membership = factor(clean_data$group_membership,
levels=c("AMC", "BSO", "BSO-E2"))
### Dunn test
PT = dunnTest(Cancer_Treatments ~ group_membership,
data=clean_data,
method="bonferroni")
PT
#MenopauseAge test
Menopause.aov<-aov(MenopauseAge~group_membership)
summary(Menopause.aov)
leveneTest(MenopauseAge~group_membership) #assess the equality of the variances
# Extract the residuals
Menopause_aovresiduals <- residuals(object = Menopause.aov )
shapiro.test(x = Menopause_aovresiduals )
#Years of Education
EduYears.aov<-aov(EduYears~group_membership)
summary(EduYears.aov)
leveneTest(EduYears~group_membership) #assess the equality of the variances
# Extract the residuals
EduYears_aovresiduals <- residuals(object = EduYears.aov )
shapiro.test(x = EduYears_aovresiduals )
#BMI
BMI.aov<-aov(BMI~group_membership)
summary(BMI.aov)
leveneTest(BMI~group_membership) #assess the equality of the variances
# Extract the residuals
BMI_aovresiduals <- residuals(object = BMI.aov )
shapiro.test(x = BMI_aovresiduals )
#P value less than 0.05 so perform Kruskal-Wallis Rank Sum Test
kruskal.test(BMI~group_membership)
#PdG
PdG.aov <- aov(PdG~group_membership)
summary(PdG.aov)
#E1G
E1G.aov<-aov(E1G~group_membership)
summary(E1G.aov)
leveneTest(E1G~group_membership) #assess the equality of the variances
# Extract the residuals
E1G_aovresiduals <- residuals(object = E1G.aov )
shapiro.test(x = E1G_aovresiduals )
#P value smaller than 0.05 so run a Kruskal-Wallis Rank Sum Test
kruskal.test(E1G~group_membership)
# Post-hoc for KW is the Dunn Test
clean_data$group_membership = factor(clean_data$group_membership,
levels=c("AMC", "BSO", "BSO-E2"))
### Dunn test
PT = dunnTest(E1G ~ group_membership,
data=clean_data,
method="bonferroni")
PT
#Make graphs
#Compute ANOVA between two independent variables, genotype and group membership
#Since the samples sizes in the groups are uneven, you ned to do an unbalanced two way ANOVA
#When you run a linear model, there are some missing values in the COMT column. Remove those missing values
clean_data$COMT2 <- factor(clean_data$COMT, exclude = "", levels = c("Met","Het","Val"))
#Remove missing values for in the BDNF column too
clean_data$BDNF2 <- factor(clean_data$BDNF, exclude = "", c("Met_carrier", "Val"))
#First, check the contrasts. What are the contrasts? It shows how one group is compared to another group.
#Each type of contrast will give you a different sum of squares. The default setting is contr.treatment(for unordered data) and contr.poly(for ordered data)
#when you open R and you check the contrasts, it'll show you treatment and poly.
#You need to set it to contr.sum and contr.poly. You can't use treatment because the tests wont give you sensible results. This is in the help file(According to the person who made it in car)
#The data is already sorted intro reference groups. The reference group is what the other groups are compared to. Use the following code to change the reference groups
clean_data$group_membership <- factor(clean_data$group_membership, levels = c("AMC","BSO-E2","BSO")) #to make AMC the reference group
clean_data$group_membership <- factor(clean_data$group_membership, levels = c("BSO","BSO-E2","AMC")) #to make BSO the reference group
knitr::kable(nice(ancova_RAVLTlearn))%>% # this code allows you to save out your table directly in a word doc
write2word(file="RAVLTlearn.doc", quiet = TRUE) # see this website, https://cran.r-project.org/web/packages/arsenal/vignettes/write2.html#introduction
# Two factorial ANOVAs -----
options()$contrasts #to check the contrasts
options(contrasts = c("contr.sum","contr.poly"))
Anova(lm(SPWM_WME_T1 ~ COMT2 * group_membership + Cancer_Treatments, data = clean_data), type = 3)
#Since results are significant, run post hoc test
#The TukeyHSD function is only for aov, it cannot be used on the Anova function
#Run a linear regression to see where the differences in the groups is
#Don't forget to change the contrasts back to treatment and poly! (the default)
options(contrasts = c("contr.treatment", "contr.poly"))
COMT_SPWM <- lm(SPWM_WME_T1 ~ COMT2 * group_membership + Cancer_Treatments, data = clean_data)
summary(COMT_SPWM)
#Results:
#COMT2Val: The difference between Val and Het is significant in AMCs
#group_membershipBSO-E2: The difference between BSO-E2 and AMCs is significant for Het
#COMT2Val:group_membershipBSO-E2: The interaction term, the difference of the differences
#My interpretation: difference between Vals in BSO and BSO-E2 (in this case the constant is the val)
#Het and AMCs are the reference groups here. To check reference groups, use contrasts function
contrasts(clean_data$group_membership)
clean_data$COMT2 <- factor(clean_data$COMT2, levels = c("Het","Met","Val"))
clean_data$group_membership <- factor(clean_data$group_membership, levels = c("AMC","BSO-E2","BSO"))
clean_data$group_membership <- factor(clean_data$group_membership, levels = c("BSO","BSO-E2","AMC"))
#Since the value for AMC is 0 in both BSO and BSO-E2, it is the reference group
contrasts(clean_data$COMT2)
#Since the value for Het is 0 in both Met and Val, it is the reference group
# Makes more sense if the reference group is Met (we want to compare the Met to Val)
#I want to know if DigitsForward scores have an effect on COMT genotype and group membership
options()$contrasts #to check the contrasts
options(contrasts = c("contr.sum","contr.poly"))
Anova(lm(DigitsForward ~ COMT2 * group_membership + Cancer_Treatments, data = clean_data), type = 3)
#Not signficant
#I want to know if DigitsBackward scores have an effect on COMT genotype and group membership
options()$contrasts #to check the contrasts
options(contrasts = c("contr.sum","contr.poly"))
Anova(lm(DigitsBackward ~ COMT2 * group_membership + Cancer_Treatments, data = clean_data), type = 3)
#Not significant
#I want to know if DigitsTotal scores have an effect on COMT genotype and group membership
options()$contrasts #to check the contrasts
options(contrasts = c("contr.sum","contr.poly"))
Anova(lm(DigitsTotal ~ COMT2 * group_membership + Cancer_Treatments, data = clean_data), type = 3)
#Not significant, but there's a trend in the COMT2 scores
#Run a linear model maybe
options(contrasts = c("contr.treatment", "contr.poly"))
COMT_DigitsTotal <- lm(DigitsTotal ~ COMT2 * group_membership + Cancer_Treatments, data = clean_data)
summary(COMT_DigitsTotal)
#I want to know if RAVLT_learn scores have an effect on BDNF genotype and group membership
#To change reference groups
clean_data$group_membership <- factor(clean_data$group_membership, levels = c("AMC","BSO-E2","BSO"))
clean_data$group_membership <- factor(clean_data$group_membership, levels = c("BSO","BSO-E2","AMC"))
#To check reference groups
contrasts(clean_data$BDNF_combined)
contrasts(clean_data$group_membership)
options()$contrasts #to check the contrasts
options(contrasts = c("contr.sum","contr.poly"))
Anova(lm(RAVLT_Learn ~ BDNF_combined * group_membership + Cancer_Treatments, data = clean_data), type = 3)
#contr.sum(clean_data$BDNF_combined)
options(contrasts = c("contr.treatment", "contr.poly"))
BDNF_RAVLT<- lm(RAVLT_Learn ~ BDNF_combined * group_membership + Cancer_Treatments, data = clean_data)
summary(BDNF_RAVLT)
#Logical Memory Immediate
options()$contrasts #to check the contrasts
options(contrasts = c("contr.sum","contr.poly"))
Anova(lm(LMA_Imm_Verbatim ~ BDNF_combined * group_membership + Cancer_Treatments, data = clean_data), type = 3)
options(contrasts = c("contr.treatment", "contr.poly"))
BDNF_LMA_Imm_RAVLT<- lm(LMA_Imm_Verbatim ~ BDNF_combined * group_membership + Cancer_Treatments, data = clean_data)
summary(BDNF_LMA_Imm_RAVLT)
#Logical Memory Delayed
options()$contrasts #to check the contrasts
options(contrasts = c("contr.sum","contr.poly"))
Anova(lm(LMA_Del_Verbatim ~ BDNF_combined * group_membership + Cancer_Treatments, data = clean_data), type = 3)
options(contrasts = c("contr.treatment", "contr.poly"))
BDNF_LMA_Del_RAVLT<- lm(LMA_Del_Verbatim ~ BDNF_combined * group_membership + Cancer_Treatments, data = clean_data)
summary(BDNF_LMA_Del_RAVLT)
#To change the reference group
clean_data$group_membership <- factor(clean_data$group_membership, levels = c("AMC","BSO-E2","BSO"))
clean_data$group_membership <- factor(clean_data$group_membership, levels = c("BSO","BSO-E2","AMC"))
#making graphs ----
#Make a graph for the E1G levels
#make box plot for E1G levels
#First make a graph theme that can be applied to all your graphs
E1GTheme <- theme(plot.title= element_text(family = "Palatino", face = "bold")) #determining font
E1GPlot <- ggplot(clean_data, aes(X=group_membership, y = E1G, fill=group_membership)) + geom_boxplot(outlier.size = .8, lwd=1)
print(E1GPlot + E1GTheme + labs(title="17-β-Estradiol Levels in Different Women Groups",
y = "Estrone-3-Glucoronide (ng/ml)", x="Group Membership",
fill = "Group Membership"))
#Box plot with jittered points
#Change outline colours by groups
#Use a custom colour palette
E1G_Boxplot <- ggboxplot(clean_data, x = "group_membership", y = "E1G",
color = "group_membership", palette =c("#00AFBB", "#E7B800", "#FC4E07"),
add = "jitter", shape = "group_membership", legend = "none")
E1GTheme <- theme(plot.title= element_text(family = "Helvetica", size = (13), hjust = 0.5))
my_comparisons <- list( c("BSO", "BSO-E2"), c("BSO-E2", "AMC"), c("BSO", "AMC") ) #if you want to include Kruskal-Wallis values
print(E1G_Boxplot
+ E1GTheme
+ stat_compare_means(comparisons = my_comparisons) #if you want to include Kruskal-Wallis values
+ stat_compare_means(label.y = 165) #if you want to include Kruskal Wallis values
+ labs(title="17-β-Estradiol Levels in Different Women Groups",
y = "Estrone-3-Glucoronide (ng/ml)", x="Group Membership",
fill = "Group Membership"))
#Violin plot
E1GViolinplot <- ggviolin(clean_data, x = "group_membership", y = "E1G", fill = "group_membership",
palette = c("#00AFBB", "#E7B800", "#FC4E07"),
add.params = list(fill = "white"),
add = "boxplot", legend = "none")
E1GTheme <- theme(plot.title= element_text(family = "Helvetica", size = (13), hjust = 0.5))
my_comparisons <- list( c("BSO", "BSO-E2"), c("BSO-E2", "AMC"), c("BSO", "AMC") ) #if you want to include Kruskal-Wallis values
print(E1GViolinplot
+ E1GTheme
+ stat_compare_means(comparisons = my_comparisons) #if you want to include Kruskal-Wallis values
+ stat_compare_means(label.y = 165) #if you want to include Kruskal Wallis values
+ labs(title="17-β-Estradiol Levels in Different Women Groups",
y = "Estrone-3-Glucoronide (ng/ml)", x="Group Membership",
fill = "Group Membership"))
#Make a bar graph for COMT and SPWM
BarTheme <- theme(plot.title= element_text(family = "Helvetica", size = (13),hjust = 0.5),
panel.background = element_rect(fill = "white", size = 4),
axis.line = element_line(size = 0.5, colour = "black"),
axis.ticks = element_line(colour = "black", size = 0.5))
SPWMBar <- ggplot(clean_data, aes(x=group_membership, y=SPWM_WME_T1, fill=COMT2)) +
stat_summary(fun=mean, geom="bar", position=position_dodge(width=.9)) +
stat_summary(fun.data=mean_se, geom="errorbar", position=position_dodge(width=.9), width=.3) + # need to set width=.9, or else will squish all the errorbars together
labs(title = "Spatial Working Memory Errors in COMT Genotypes", x = "Group Membership", y = "Working Memory Errors") +
scale_fill_discrete(name = "COMT Genotype", labels = c("Val/Met", "Met/Met", "Val/Val"))
print(SPWMBar + BarTheme)
#?ggpar
#Make a bar graph for COMT and DigitsForward
DigitsForwardBar <- ggplot(clean_data, aes(x=group_membership, y=DigitsForward, fill=COMT2)) +
stat_summary(fun=mean, geom="bar", position=position_dodge(width=.9)) +
stat_summary(fun.data=mean_se, geom="errorbar", position=position_dodge(width=.9), width=.3) + # need to set width=.9, or else will squish all the errorbars together
labs( x = "Group Membership", y = "Mean Score") +
scale_fill_discrete(name = "COMT Genotype", labels = c("Val/Met", "Met/Met", "Val/Val"))
print(DigitsForwardBar + BarTheme)
#Make a bar graph for COMT and DigitsBackwards
DigitsBackwardBar <- ggplot(clean_data, aes(x=group_membership, y=DigitsBackward, fill=COMT2)) +
stat_summary(fun=mean, geom="bar", position=position_dodge(width=.9)) +
stat_summary(fun.data=mean_se, geom="errorbar", position=position_dodge(width=.9), width=.3) + # need to set width=.9, or else will squish all the errorbars together
labs(x = "Group Membership", y = "Mean Score") +
scale_fill_discrete(name = "COMT Genotype", labels = c("Val/Met", "Met/Met", "Val/Val"))
print(DigitsBackwardBar + BarTheme)
#Make a bar graph for BDNF and RAVLT_Learn
RAVLT_LearnBar <- ggplot(clean_data, aes(x=group_membership, y=RAVLT_Learn, fill=BDNF_combined)) +
stat_summary(fun=mean, geom="bar", position=position_dodge(width=.9)) +
stat_summary(fun.data=mean_se, geom="errorbar", position=position_dodge(width=.9), width=.3) + # need to set width=.9, or else will squish all the errorbars together
labs(x = "Group Membership", y = "RAVLT Percent Learn (A5/A6 x 100%)") +
scale_fill_discrete(name = "BDNF Genotype", labels = c("Met Carrier", "Val/Val"))
print(RAVLT_LearnBar + BarTheme)
#Make a bar graph for BDNF and LM_Imm_Verbatim
LMA_Imm_VerbatimBar <- ggplot(clean_data, aes(x=group_membership, y=LMA_Imm_Verbatim, fill=BDNF_combined)) +
stat_summary(fun=mean, geom="bar", position=position_dodge(width=.9)) +
stat_summary(fun.data=mean_se, geom="errorbar", position=position_dodge(width=.9), width=.3) + # need to set width=.9, or else will squish all the errorbars together
labs( x = "Group Membership", y = "Immediate Verbatim Mean Score") +
scale_fill_discrete(name = "BDNF Genotype", labels = c("Met Carrier", "Val/Val"))
print(LMA_Imm_VerbatimBar + BarTheme)
#Make a bar graph for BDNF and LM_Imm_Verbatim
LMA_Del_VerbatimBar <- ggplot(clean_data, aes(x=group_membership, y=LMA_Del_Verbatim, fill=BDNF_combined)) +
stat_summary(fun=mean, geom="bar", position=position_dodge(width=.9)) +
stat_summary(fun.data=mean_se, geom="errorbar", position=position_dodge(width=.9), width=.3) + # need to set width=.9, or else will squish all the errorbars together
labs(x = "Group Membership", y = "Delayed Verbatim Mean Score") +
scale_fill_discrete(name = "BDNF Genotype", labels = c("Met Carrier", "Val/Val"))
print(LMA_Del_VerbatimBar + BarTheme)
#Hardy Weinberg Test for COMT AMCs
x <- c(MM=4,VM=11,VV=7)
HW.test <- HWChisq(x, verbose=TRUE)
#P value 0.79 therefore population in HWE, X2 = 0.0689
#Hardy Weinberg Test for COMT BSOs
x <- c(MM=6,VM=6,VV=5)
HW.test <- HWChisq(x, verbose=TRUE)
#P value 0.39 therefore population in HWE, X2 = 0.7155
#Hardy Weinberg Test for COMT BSO-E2s
x <- c(MM=4,VM=7,VV=4)
HW.test <- HWChisq(x, verbose=TRUE)
#P value 0.85 therefore population in HWE, X2 = 0.03
#Hardy Weinberg Test for BDNF AMCs
x <- c(MM=1,VM=3,VV=9)
HW.test <- HWChisq(x, verbose=TRUE)
#P value 0.78 therefore population in HWE, X2 = 0.07
#Hardy Weinberg Test for BDNF BSOs
x <- c(MM=2,VM=5,VV=15)
HW.test <- HWChisq(x, verbose=TRUE)
#P value 0.3 therefore population in HWE, X2 = 0.77
#Hardy Weinberg Test for BDNF BSOE2s
x <- c(MM=0,VM=2,VV=11)
HW.test <- HWChisq(x, verbose=TRUE)
#P value is 0.12 therefore population in HWE, X2 = 2.4
help(HardyWeinberg)
citation("car")
citation("HardyWeinberg")