-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathcode_based_on_simulated_data.r
148 lines (134 loc) · 3.24 KB
/
code_based_on_simulated_data.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
#############################
setwd(dirname(file.choose()))
read.csv("coe.csv")->coe
#count zero questions
str(coe)
size<-20
count1<-rep(-1,size)
count2<-rep(-1,size)
count3<-rep(-1,size)
coe[c(2:81),]->try1
coe[c(82:161),]->try2
coe[c(162:241),]->try3
for (i in 2:(size+1)){
sum(try1[,i]==0)/2->count1[i-1]
}
for (i in 2:(size+1)){
sum(try2[,i]==0)/2->count2[i-1]
}
for (i in 2:(size+1)){
sum(try3[,i]==0)/2->count3[i-1]
}
summary(count1)
summary(count2)
summary(count3)
count1
count2
count3
sd(count1)->sd1
mean(count1)+1.96*sd1
mean(count1)-1.96*sd1
sd(count2)->sd2
mean(count2)+1.96*sd2
mean(count2)-1.96*sd2
sd(count3)->sd3
mean(count3)+1.96*sd3
mean(count3)-1.96*sd3
############### Nov 10, compare to theoretic prob.
M<-20
AUC<-NULL
#nlam is the number of values of tuning variable
nlam<-148
for (j in 1:M){
read.csv(paste('pre_link',j,'.csv',sep=''))->pre
read.csv(paste('sim',j,'_da',1,'.csv',sep=''))->da1
read.csv(paste('sim',j,'_da',2,'.csv',sep=''))->da2
read.csv(paste('sim',j,'_da',3,'.csv',sep=''))->da3
y<-da1$y
n<-ncol(da3)
r<-ncol(da2)
row<-nrow(da2)
nfarm<-row-1
AUC1<-rep(0,nlam)
#-------------------------------------------
for (i in 1:148){
AUC1[i]<-rocTest(y,x=list(Model1=y,Model0=pre[,i]),L=matrix(c(1,-1),1))$th[2]
}
AUC<-rbind(AUC,AUC1)
}
write.csv(AUC,"AUC.csv")
############### Nov 11
read.csv("AUC.csv")->auc
auc[,c(2:149)]->auc
str(auc)
r<-nrow(auc)
idx<-rep(0,r)
maxauc<-rep(0,r)
for (i in 1:r){
maxauc[i]<-max(auc[i,][auc[i,]==max(auc[i,])])
}
###############calculate the theoretical value of prob.#####
M<-20
link<-NULL
theory_auc<-rep(0,M)
for (j in 1:M){
read.csv(paste('dummy.sim',j,'.csv',sep=''))->dummy.sim1
read.csv(paste('sim',j,'_da',1,'.csv',sep=''))->da1
y<-da1$y
c(rep(c(1/4,0,-1/4),40),rep(c(1/4,0,0),40),rep(c(0,0,0),40))->s1
as.matrix(dummy.sim1)%*%s1-40/3/4->link1
theory_auc[j]<-rocTest(y,x=list(Model1=y,Model0=link1),L=matrix(c(1,-1),1))$th[2]
link<-cbind(link,link1)
}
data.frame(link)->link
data.frame(rbind(theory_auc,maxauc))->compare_auc
compare_auc
#str(compare_auc)
write.csv(link,'the_link.csv')
write.csv(compare_auc,"compare_auc.csv")
##################
############### Feb 21,2012
M<-20
AUC<-NULL
#read.csv('compare_auc.csv',row.names=1)->auc
#nlam is the number of values of tuning variable
for (j in 1:M){
read.csv(paste('sim',j,'_da1_SAS','.csv',sep=''))->pre
y<-pre$y
#read.csv(paste('sim',j,'_da',1,'.csv',sep=''))->da1
#read.csv(paste('sim',j,'_da',2,'.csv',sep=''))->da2
#read.csv(paste('sim',j,'_da',3,'.csv',sep=''))->da3
AUC[j]<-rocTest(y,x=list(Model1=y,Model0=pre$phat),L=matrix(c(1,-1),1))$th[2]
}
str(auc)
rbind(auc,sas_auc=AUC)->auc
t.test(auc[2,],auc[3,])
write.csv(t(AUC),"AUC.csv")
write.csv(auc,"compare_auc.csv")
############### Nov 11
read.csv("AUC.csv")->auc
auc[,c(2:149)]->auc
str(auc)
r<-nrow(auc)
idx<-rep(0,r)
maxauc<-rep(0,r)
for (i in 1:r){
maxauc[i]<-max(auc[i,][auc[i,]==max(auc[i,])])
}
#####wilcox rank test
read.csv(file.choose())->auc
auc2<-auc[,c(2:21)]
x0.1<-auc2[1,]
y0.1<-auc2[2,]
x0.25<-auc2[3,]
y0.25<-auc2[4,]
x0.5<-auc2[5,]
y0.5<-auc2[6,]
x1<-auc2[7,]
y1<-auc2[8,]
x2<-auc2[9,]
y2<-auc2[10,]
mean(t(y0.1))
sd(t(y0.1))
str(as.numeric(x0.1))
wilcox.test(as.numeric(y2) - as.numeric(x2))