-
Notifications
You must be signed in to change notification settings - Fork 7
/
Copy path.Rhistory
242 lines (242 loc) · 9.38 KB
/
.Rhistory
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
devtools::document('scLink')
setwd("~/Dropbox/Rpkgs-dev/sclink_dev/scLink")
devtools::document('scLink')
setwd("~/Dropbox/Rpkgs-dev/sclink_dev/")
devtools::document('scLink')
system('R CMD Rd2pdf scLink')
devtools::document('scLink')
devtools::document('scLink')
system('R CMD Rd2pdf scLink')
rm(list = ls())
install.packages("~/Dropbox/Rpkgs-dev/sclink_dev/scLink_0.0.1.tar.gz", repos = NULL)
library(scLink)
count = system.file("extdata", "example.rds", package = "scLink")
count = system.file("extdata", "example.rds", package = "scLink")
genes = system.file("extdata", "genes.rds", package = "scLink")
count.norm = sclink_norm(count, scale.factor = 1e6, filter.genes = FALSE, gene.names = genes)
count[1:3,1:3]
count = readRDS(system.file("extdata", "example.rds", package = "scLink"))
genes = readRDS(system.file("extdata", "genes.rds", package = "scLink"))
count.norm = sclink_norm(count, scale.factor = 1e6, filter.genes = FALSE, gene.names = genes)
corr = sclink_cor(expr = count.norm, ncores = 1)
corr[1:3,1:3]
plot(count[,"Rn45s"], count[,"Eef1a1"])
plot(count.norm[,"Rn45s"], count.norm[,"Eef1a1"])
networks = sclink_net(expr = count.norm, ncores = 1, lda = seq(1, 0.1, -0.05))
networks$cor[1:3,1:3]
net1 = networks$summary[[1]]
net1$nedge
net1$lambda
sapply(networks$summary, function(x) x$lambda)
sapply(networks$summary, function(x) x$nedge)
sapply(networks$summary, function(x) x$bic)
rm(list = ls())
#install.packages("~/Dropbox/Rpkgs-dev/sclink_dev/scLink_0.0.1.tar.gz", repos = NULL)
library(devtools)
install_github("Vivianstats/scLink")
library(scLink)
count = readRDS(system.file("extdata", "example.rds", package = "scLink"))
genes = readRDS(system.file("extdata", "genes.rds", package = "scLink"))
count.norm = sclink_norm(count, scale.factor = 1e6, filter.genes = FALSE, gene.names = genes)
corr = sclink_cor(expr = count.norm, ncores = 1)
corr[1:3,1:3]
networks = sclink_net(expr = count.norm, ncores = 1, lda = seq(1, 0.1, -0.05))
sapply(networks$summary, function(x) x$lambda)
sapply(networks$summary, function(x) x$nedge)
sapply(networks$summary, function(x) x$bic)
240*0.75/0.91
200/240
install.packages("pwr")
library(pwr)
pwr.t2n.test(n1 = 205, n2= 463, d = 1.6, sig.level = 0.05,
alternative = "two.sided")
install.packages("powerAnalysis")
library(powerAnalysis)
ES.t.two(m1 = 110.1, m2 = 111.7, sd1 = 23.8, sd2 = 25.8,
n1 = 205, n2 = 463, alternative = "two.sided")
pwr.t2n.test(n1 = 205, n2= 463, d = 0.063, sig.level = 0.05,
alternative = "two.sided")
pwr.t2n.test(n1 = c(205, 210), n2= c(463,480), d = 0.063, sig.level = 0.05,
alternative = "two.sided")
setwd("~/Box/Collaboration/201909Adana/sample-size")
prior = readxlsx("IHC-summary.xlsx")
library(openxlsx)
prior = readxlsx("IHC-summary.xlsx")
prior = read.xlsx("IHC-summary.xlsx")
prior
i = 1
es = ES.t.two(m1 = prior[i,"mu1"], m2 = prior[i,"mu2"],
sd1 = prior[i,"sd1"], sd2 = prior[i,"sd2"],
n1 = prior[i,"n1"], n2 = prior[i,"n2"], alternative = "two.sided")
es$d
NN = seq(500,2000,100)
NN
NN = seq(500,2000,100)
n1 = round(NN * prior[i,"mu1"] / (prior[i,"mu1"] + prior[i,"mu2"]))
n2 = NN - n1
test = pwr.t2n.test(n1 = n1, n2 = n2, d = es$d, sig.level = 0.05,
alternative = "two.sided")
test
test$power
plot(NN, test$power)
plot(NN, test$power)
plot(NN, test$power, type = "b")
plot(NN, test$power, type = "b")
res = lapply(1:nrow(prior), function(i){
es = ES.t.two(m1 = prior[i,"mu1"], m2 = prior[i,"mu2"],
sd1 = prior[i,"sd1"], sd2 = prior[i,"sd2"],
n1 = prior[i,"n1"], n2 = prior[i,"n2"], alternative = "two.sided")
n1 = round(NN * prior[i,"mu1"] / (prior[i,"mu1"] + prior[i,"mu2"]))
n2 = NN - n1
test = pwr.t2n.test(n1 = n1, n2 = n2, d = es$d, sig.level = 0.05,
alternative = "two.sided")
data.frame(size = NN, power = test$power, marker = prior[i,"marker"],
status1 = prior[i,"status1"], status2 = prior[i,"status2"])
})
res = Reduce(rbind, res)
library(ggplot2); theme_set(theme_classic())
ggplot(res, aes(x = size, y = power))+
geom_point() + geom_line() +
facet_grid(~marker)
NN = seq(500,2000,100)
res = lapply(1:nrow(prior), function(i){
es = ES.t.two(m1 = prior[i,"mu1"], m2 = prior[i,"mu2"],
sd1 = prior[i,"sd1"], sd2 = prior[i,"sd2"],
n1 = prior[i,"n1"], n2 = prior[i,"n2"], alternative = "two.sided")
n1 = round(NN * 0.5)
n2 = NN - n1
test = pwr.t2n.test(n1 = n1, n2 = n2, d = es$d, sig.level = 0.05,
alternative = "two.sided")
data.frame(size = NN, power = test$power, marker = prior[i,"marker"],
status1 = prior[i,"status1"], status2 = prior[i,"status2"])
})
res = Reduce(rbind, res)
ggplot(res, aes(x = size, y = power))+
geom_point() + geom_line() +
facet_grid(~marker)
NN = seq(600,1800,100)
res = lapply(1:nrow(prior), function(i){
es = ES.t.two(m1 = prior[i,"mu1"], m2 = prior[i,"mu2"],
sd1 = prior[i,"sd1"], sd2 = prior[i,"sd2"],
n1 = prior[i,"n1"], n2 = prior[i,"n2"], alternative = "two.sided")
n1 = round(NN * prior[i,"mu1"] / (prior[i,"mu1"] + prior[i,"mu2"]))
n2 = NN - n1
test = pwr.t2n.test(n1 = n1, n2 = n2, d = es$d, sig.level = 0.05,
alternative = "two.sided")
data.frame(size = NN, power = test$power, marker = prior[i,"marker"],
status1 = prior[i,"status1"], status2 = prior[i,"status2"])
})
res = Reduce(rbind, res)
ggplot(res, aes(x = size, y = power))+
geom_point() + geom_line() +
facet_grid(~marker)
i = 4
es = ES.t.two(m1 = prior[i,"mu1"], m2 = prior[i,"mu2"],
sd1 = prior[i,"sd1"], sd2 = prior[i,"sd2"],
n1 = prior[i,"n1"], n2 = prior[i,"n2"], alternative = "two.sided")
es
prior[i,"mu1"]
prior[i,"mu2"]
prior = read.xlsx("IHC-summary.xlsx")
NN = seq(600,1800,100)
res = lapply(1:nrow(prior), function(i){
es = ES.t.two(m1 = prior[i,"mu1"], m2 = prior[i,"mu2"],
sd1 = prior[i,"sd1"], sd2 = prior[i,"sd2"],
n1 = prior[i,"n1"], n2 = prior[i,"n2"], alternative = "two.sided")
n1 = round(NN * prior[i,"mu1"] / (prior[i,"mu1"] + prior[i,"mu2"]))
n2 = NN - n1
test = pwr.t2n.test(n1 = n1, n2 = n2, d = es$d, sig.level = 0.05,
alternative = "two.sided")
data.frame(size = NN, power = test$power, marker = prior[i,"marker"],
status1 = prior[i,"status1"], status2 = prior[i,"status2"])
})
res = Reduce(rbind, res)
ggplot(res, aes(x = size, y = power))+
geom_point() + geom_line() +
facet_grid(~marker)
NN = seq(600,6000,500)
res = lapply(1:nrow(prior), function(i){
es = ES.t.two(m1 = prior[i,"mu1"], m2 = prior[i,"mu2"],
sd1 = prior[i,"sd1"], sd2 = prior[i,"sd2"],
n1 = prior[i,"n1"], n2 = prior[i,"n2"], alternative = "two.sided")
n1 = round(NN * prior[i,"mu1"] / (prior[i,"mu1"] + prior[i,"mu2"]))
n2 = NN - n1
test = pwr.t2n.test(n1 = n1, n2 = n2, d = es$d, sig.level = 0.05,
alternative = "two.sided")
data.frame(size = NN, power = test$power, marker = prior[i,"marker"],
status1 = prior[i,"status1"], status2 = prior[i,"status2"])
})
res = Reduce(rbind, res)
ggplot(res, aes(x = size, y = power))+
geom_point() + geom_line() +
facet_grid(~marker)
NN = seq(600,1800,200)
res = lapply(1:nrow(prior), function(i){
es = ES.t.two(m1 = prior[i,"mu1"], m2 = prior[i,"mu2"],
sd1 = prior[i,"sd1"], sd2 = prior[i,"sd2"],
n1 = prior[i,"n1"], n2 = prior[i,"n2"], alternative = "two.sided")
n1 = round(NN * prior[i,"mu1"] / (prior[i,"mu1"] + prior[i,"mu2"]))
n2 = NN - n1
test = pwr.t2n.test(n1 = n1, n2 = n2, d = es$d, sig.level = 0.05,
alternative = "two.sided")
data.frame(size = NN, power = test$power, marker = prior[i,"marker"],
status1 = prior[i,"status1"], status2 = prior[i,"status2"])
})
res = Reduce(rbind, res)
ggplot(res, aes(x = size, y = power))+
geom_point() + geom_line() +
facet_grid(~marker)
paste0(res$status1, res$status2)
paste0(res$status1, res$status2, collapse = " vs. ")
paste0(res$status1, res$status2, sep = " vs. ")
res$study = paste0(res$status1, sep = " vs. ", res$status2)
ggplot(res, aes(x = size, y = power))+
geom_point() + geom_line() +
facet_grid(study~marker)
prior = read.xlsx("IHC-summary.xlsx")
NN = seq(600,1800,200)
res = lapply(1:nrow(prior), function(i){
es = ES.t.two(m1 = prior[i,"mu1"], m2 = prior[i,"mu2"],
sd1 = prior[i,"sd1"], sd2 = prior[i,"sd2"],
n1 = prior[i,"n1"], n2 = prior[i,"n2"], alternative = "two.sided")
n1 = round(NN * prior[i,"mu1"] / (prior[i,"mu1"] + prior[i,"mu2"]))
n2 = NN - n1
test = pwr.t2n.test(n1 = n1, n2 = n2, d = es$d, sig.level = 0.05,
alternative = "two.sided")
data.frame(size = NN, power = test$power, marker = prior[i,"marker"],
status1 = prior[i,"status1"], status2 = prior[i,"status2"])
})
res = Reduce(rbind, res)
res$study = paste0(res$status1, sep = " vs. ", res$status2)
ggplot(res, aes(x = size, y = power))+
geom_point() + geom_line() +
facet_grid(study~marker)
ggplot(res, aes(x = size, y = power))+
geom_point() + geom_line() +
facet_wrap(study~marker)
ggplot(res, aes(x = size, y = power))+
geom_point() + geom_line() +
facet_wrap(study~marker, ncol = 5)
library(ggplot2); theme_set(theme_bw())
ggplot(res, aes(x = size, y = power))+
geom_point() + geom_line() +
facet_grid(study~marker)
prior = read.xlsx("IHC-summary.xlsx")
NN = seq(600,1800,200)
res = lapply(1:nrow(prior), function(i){
es = ES.t.two(m1 = prior[i,"mu1"], m2 = prior[i,"mu2"],
sd1 = prior[i,"sd1"], sd2 = prior[i,"sd2"],
n1 = prior[i,"n1"], n2 = prior[i,"n2"], alternative = "two.sided")
n1 = round(NN * prior[i,"mu1"] / (prior[i,"mu1"] + prior[i,"mu2"]))
n2 = NN - n1
test = pwr.t2n.test(n1 = n1, n2 = n2, d = es$d, sig.level = 0.05,
alternative = "two.sided")
data.frame(size = NN, power = test$power, marker = prior[i,"marker"],
status1 = prior[i,"status1"], status2 = prior[i,"status2"])
})
res = Reduce(rbind, res)
res$study = paste0(res$status1, sep = " vs. ", res$status2)
ggplot(res, aes(x = size, y = power))+
geom_point() + geom_line() +
facet_grid(study~marker)
ggsave("samplesize-0305.pdf", width = 8, height = 6)