-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathgeoMatch_simTests.R
87 lines (75 loc) · 4.73 KB
/
geoMatch_simTests.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
#SciClone Edition
#library(devtools)
#install_github('itpir/geoMatch')
library(geoMatch)
source("/sciclone/home00/geogdan/geoMatch_testing/SimTests/spatial.simulation.R")
#For SciClone - LEAVE THIS AT 1! Iterations are handled in python.
iterations <- 1
record.dataframe <- data.frame(n.points=round(runif(iterations,min=200,max=200)),
cov.Decay = runif(iterations, min=0.01, max=0.01),
se.Decay = runif(iterations, min=0.01, max=0.01),
t.Decay = runif(iterations, min=50, max=50),
sim.T.e = runif(iterations, min=0.0, max=0.0),
T.percent = runif(iterations, min=0.5, max=0.5),
sim.Y.scale = runif(iterations, min=0.0, max=0.0),
Theta = runif(iterations, min=1.0, max=1.0),
sim.Y.cov.effect = runif(iterations, min=0.0, max=0.0),
sim.Y.het.effect = runif(iterations, min=0.0, max=0.0),
sim.Y.e = runif(iterations, min=0.0, max=0.0),
spill.mag = runif(iterations, min=0.0, max=2.0),
avg.spill = NA,
est.T.lm = NA,
est.T.matchit = NA,
est.T.geoMatch = NA,
abs.error.lm = NA,
abs.error.matchit = NA,
abs.error.geoMatch = NA,
mean.error.lm = NA,
mean.error.matchit = NA,
mean.error.geoMatch = NA,
true.T.ie = NA,
converge = NA)
for(i in 1:iterations)
{
print(i)
it.spdf <- spatial.simulation(n.points = record.dataframe["n.points"][i,],
cov.Decay = record.dataframe["cov.Decay"][i,],
se.Decay = record.dataframe["se.Decay"][i,],
t.Decay = record.dataframe["t.Decay"][i,],
sim.T.e= record.dataframe["sim.T.e"][i,],
T.percent = record.dataframe["T.percent"][i,],
sim.Y.scale = record.dataframe["sim.Y.scale"][i,],
Theta = record.dataframe["Theta"][i,],
sim.Y.cov.effect = record.dataframe["sim.Y.cov.effect"][i,],
sim.Y.het.effect = record.dataframe["sim.Y.het.effect"][i,],
sim.Y.e = record.dataframe["sim.Y.e"][i,],
spill.mag = record.dataframe["spill.mag"][i,])
it.geoMatch <- geoMatch(T ~ X,
method = "nearest",
caliper=0.25,
data = it.spdf,
outcome.variable="Y",
outcome.suffix="_adjusted",
max.it=10,
report = TRUE)
record.dataframe$converge <- it.geoMatch$optim
it.Match <- matchit(T ~ X,
method = "nearest",
caliper=0.25,
data = it.spdf@data)
lm.model <- lm(Y ~ T + X, data=it.spdf@data)
matchit.model <- lm(Y ~ T + X, data=match.data(it.Match))
geo.model <- lm(Y_adjusted ~ T + X, data=match.data(it.geoMatch))
record.dataframe["mean.error.lm"][i,] <- mean((lm.model$coefficients["T"] - mean(it.spdf@data[it.spdf@data$T == 1,]["ie.spill"][[1]])))
record.dataframe["mean.error.matchit"][i,] <- mean((matchit.model$coefficients["T"] - mean(it.spdf@data[it.spdf@data$T == 1,]["ie.spill"][[1]])))
record.dataframe["mean.error.geoMatch"][i,] <- mean((geo.model$coefficients["T"] - mean(it.spdf@data[it.spdf@data$T == 1,]["ie.spill"][[1]])))
record.dataframe["abs.error.lm"][i,] <- mean(abs((lm.model$coefficients["T"] - it.spdf@data[it.spdf@data$T == 1,]["ie.spill"][[1]])))
record.dataframe["abs.error.matchit"][i,] <- mean(abs((matchit.model$coefficients["T"] - it.spdf@data[it.spdf@data$T == 1,]["ie.spill"][[1]])))
record.dataframe["abs.error.geoMatch"][i,] <- mean(abs((geo.model$coefficients["T"] - it.spdf@data[it.spdf@data$T == 1,]["ie.spill"][[1]])))
record.dataframe["est.T.lm"][i,] <- lm.model$coefficients["T"]
record.dataframe["est.T.matchit"][i,] <- matchit.model$coefficients["T"]
record.dataframe["est.T.geoMatch"][i,] <- geo.model$coefficients["T"]
record.dataframe["true.T.ie"][i,] <- mean(it.spdf@data[it.spdf@data$T == 1,]["ie.spill"][[1]])
record.dataframe["avg.spill"][i,] <- mean(it.spdf@data[it.spdf@data$T == 0,]["t.spill"][[1]])
}
write.csv(record.dataframe, paste('/sciclone/home00/geogdan/geoMatch_testing/Results/out_',date(),runif(1),'.csv',sep=""))