/
example.R
57 lines (41 loc) · 1.65 KB
/
example.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
# Author: Yen-Chi Chen, Christopher R. Genovese, Larry Wasserman
# Maintainer: Yen-Chi Chen <ga014528@gmail.com>
# Reference: Statistical Inference using Morse-Smale Complex
# Date: 06/22/2015
source("MSHD.R")
library(mclust) #for GvHD dataset
####
data(GvHD)
D0 = GvHD.control
D1 = GvHD.pos
#### parameters
h = 50 # smoothing
n_g = 21 # grid number
####
D_pull =rbind(D0,D1)
#### grid points
s1 = seq(from=quantile(D_pull[,1], 0.1), to=quantile(D_pull[,1], 0.9), length.out=n_g)
s2 = seq(from=quantile(D_pull[,2], 0.1), to=quantile(D_pull[,2], 0.9), length.out=n_g)
s3 = seq(from=quantile(D_pull[,3], 0.1), to=quantile(D_pull[,3], 0.9), length.out=n_g)
s4 = seq(from=quantile(D_pull[,4], 0.1), to=quantile(D_pull[,4], 0.9), length.out=n_g)
G0 = expand.grid(s1,s2,s3,s4)
#### density estimate
G0_den = kde(D0, G0, h)
G1_den = kde(D1, G0, h)
G_diff = G0_den-G1_den
G_diff_s = G_diff/(max(max(G_diff), abs(min(G_diff))))
#### High dimensional visualization
GvHD_MSHD = MSHD(G0, G_diff_s)
### visualization using modes and minima (signatures)
plot.tmp = plot(GvHD_MSHD)
text(plot.tmp$modes.vis, labels = 1:nrow(plot.tmp$modes.vis), cex=3)
text(plot.tmp$minima.vis, labels = 1:nrow(plot.tmp$minima.vis), cex=2)
### visualization using cells
plot.tmp = cell.plot(GvHD_MSHD, r_circle=1, col.cell.filled = "limegreen")
text(plot.tmp$cell.vis, labels = 1:nrow(plot.tmp$cell.vis), cex=2)
### visualizing the significant regions
a = quantile(G_diff_s, 0.9)
plot.tmp = cell.plot(GvHD_MSHD, lv.pos=a, lv.neg = -a, r_circle=1, col.cell.filled = "white")
text(plot.tmp$cell.vis, labels = 1:nrow(plot.tmp$cell.vis), cex=2, pos=2)
### energy test
E.test = MSEtest(D0,D1,G0)