dushoff.github.io/kernel.R at master · dushoff/dushoff.github.io · GitHub

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

# g, b are data frames with time and density
# L is a data frame with time, strength and something about hazard

b0fun <- function(g, r){
	step <- with(g, time[[2]]-time[[1]])
	d0 <- with(g, density*exp(-r*time))
	return(with(g, data.frame(time
		, density=d0/sum(step*d0)
	)))
}

# How much does a given r _reduce_ R for distribution g?
r2R <- function(g, r){
	step <- with(g, time[[2]]-time[[1]])
	return(with(g, 1/(step*sum(density*exp(-r*time)))))
}

R2r <- function(g, R){
	step <- with(g, time[[2]]-time[[1]])
	return(uniroot(
		function(r, g){r2R(g, r)-R}
		, interval=c(0, 1/step)
		, g=g
	)$root)
}

# Calculate the arithmetic mean inside the implicit equation for phi
phiBar <- function(b, L, phi){
	step <- with(g, time[[2]]-time[[1]])
	return(with(cbind(b,L), {
		step*sum(density*exp(phi*time)/strength)
	}))
}

phiFun <- function(b, L){
	step <- with(g, time[[2]]-time[[1]])
	return(uniroot(
		function(phi, b, L){phiBar(b, L, phi)-1}
		, interval=c(0, 1/step)
		, b=b, L=L
	)$root)
}

theFun <- function(g, L){
	step <- with(g, time[[2]]-time[[1]])
	return(with(cbind(g, L), {
		1/(step*sum(density/strength))
	}))
}

densPlot <- function(time, k, scen, color="black"){
	with(c(scen), {
		plot(time, k
			, main = disease
			, xlab = paste0("time (", unit, ")")
			, ylab = paste0("Density (1/", unit, ")")
			, type = "l", col=color
		)
	})
}

Rrplot <- function(g, R, scen){
	r <- R2r(g, R)
	b0 <- b0fun(g, r)
	print(b0)
	with(g, {
		densPlot(time, R*density, scen)

		densPlot(time, R*density, scen)
		lines(time, density, col="blue")
		lines(time, b0$density, col="red")

		densPlot(time, R*density, scen)
		lines(time, density, col="blue")

		densPlot(time, R*density, scen)
		lines(time, b0$density, col="red")
	})
}

dualPlot <- function(time, interv, intname, gen, stat, dtype){
	op <- par(mfrow=c(2, 1)
		, oma = c(0,0,0,0) + 0.1
		, mar = c(0,4,1,1) + 0.1
	)
	plot(time, interv, type="l"
		, axes=FALSE, xlab="", ylab=intname
	)
	axis(2)
	abline(h=stat, col="blue")
	par(mar = c(5,4,1,1) + 0.1)
	with(gen, plot(time, density
		, type="l", main=dtype, xlab="time (months)"
	))
	par(op)
}

intPlot <- function(time, density, strength, R){
	op <- par(mfrow=c(2, 1)
		, oma = c(0,0,0,0) + 0.1
		, mar = c(0,4,1,1) + 0.1
	)
	plot(time, strength, type="l"
		, axes=FALSE, xlab="", main="Intervention"
	)
	axis(2)

	par(mar = c(5,4,1,1) + 0.1)
	plot(time, R*density, ylab="kernel", type="l", main="Before and after")
	lines(time, R*density/strength, type="l")
	par(op)
}