forked from amymariemason/nlmr
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path3c_test_sqrt.R
72 lines (57 loc) · 1.95 KB
/
3c_test_sqrt.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
###
#Author: Amy Mason
# Purpose: test quadratic simulated data
# Date: Jan 2019
##
setwd("C:/Users/Amy/Documents/Non-linear MR/Matt Arnold Code")
require(MASS)
require(methods)
require(parallel)
require(metafor)
require(ggplot2)
require(matrixStats)
require(survival)
require(ggplot2)
require(gridExtra)
source("nlmr_functions.r")
source("nlme_summ_aes MA.r")
# fixing random numbers for repetition of what is generated
set.seed(4743045)
# creating random underlying data
###########################################################
# test nlme_summ_aes on summerised data (1 set)
beta_set1<-1.5
beta_set2<-2
lotsofdata<- create_summary_data(Ytype = "sqrt", keep = TRUE, N = 10000,
beta1 = beta_set1, beta2 = beta_set2,
quantiles = 100, confound = 0)
testdata<-lotsofdata$summary
alldata<-lotsofdata$alldata
# plot the underlying data entire
# by genetic type
ggplot(data = alldata, aes(x = X, y = "sqrt.Y") )+
geom_jitter(alpha = 0.3, aes(colour = as.factor(g)))
# underlying data check
alldata$xsqrt<-sqrt(alldata$X)
lm(data = alldata, sqrt.Y~X+xsqrt)
# non summerised fracpoly
keep1 <- fracpoly_mr(alldata$sqrt.Y, alldata$X, alldata$g, family = "gaussian",
q = 10, d = 1, fig = T)
# summerised
keep2 <- frac_poly_summ_mr(bx = testdata$BetaXG, bxse = testdata$seBetaXG,
by = testdata$BetaYG, byse = testdata$seBetaYG,
xmean = testdata$meanX, family ="gaussian",
fig = TRUE, d = 1)
summary.frac_poly_mr(keep2)
f <- function(x) (beta_set1*sqrt(x))
plot1 <- keep1$figure+ stat_function(fun = f, colour = "green")+
labs(title="fracpoly_mr")
plot2 <- keep2$figure+ stat_function(fun = f, colour = "green") +
labs(title="fracpoly_summ_mr")
grid.arrange(plot1, plot2, ncol=2)
#######
# linear.Y = beta1*X
# quad Y = <-beta2*X^2+beta1*X
# sqrt.Y = beta1*sqrt(X)
# log.Y = beta1*log(X)
# threshold.Y<-ifelse(X>beta2,beta1*X,0)