demonstration of monotonic smoothing with 2-D integral IVs

monotonic_smoothing.R

# demonstration of monotonic smoothing with 2-D integral IVs
# Harlan Harris
# harlan@harris.name

library(ggplot2)
library(locfit)
library(monoProc)

df <-
structure(list(counts = c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 
10, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 0, 1, 2, 3, 4, 5, 6, 7, 
8, 9, 10, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 0, 1, 2, 3, 4, 5, 
6, 7, 8, 9, 10, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 0, 1, 2, 3, 
4, 5, 6, 7, 8, 9, 10, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 0, 1, 
2, 3, 4, 5, 6, 7, 8, 9, 10, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 
0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 0, 1, 2, 3, 4, 5, 6, 7, 8, 
9, 10, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 0, 1, 2, 3, 4, 5, 6, 
7, 8, 9, 10, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 0, 1, 2, 3, 4, 
5, 6, 7, 8, 9, 10, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 0, 1, 2, 
3, 4, 5, 6, 7, 8, 9, 10, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 0, 
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 
10, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 0, 1, 2, 3, 4, 5, 6, 7, 
8, 9, 10, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 0, 1, 2, 3, 4, 5, 
6, 7, 8, 9, 10, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 0, 1, 2, 3, 
4, 5, 6, 7, 8, 9, 10, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 0, 1, 
2, 3, 4, 5, 6, 7, 8, 9, 10, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 
0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10), time = c(-30, -30, -30, -30, 
-30, -30, -30, -30, -30, -30, -30, -29, -29, -29, -29, -29, -29, 
-29, -29, -29, -29, -29, -28, -28, -28, -28, -28, -28, -28, -28, 
-28, -28, -28, -27, -27, -27, -27, -27, -27, -27, -27, -27, -27, 
-27, -26, -26, -26, -26, -26, -26, -26, -26, -26, -26, -26, -25, 
-25, -25, -25, -25, -25, -25, -25, -25, -25, -25, -24, -24, -24, 
-24, -24, -24, -24, -24, -24, -24, -24, -23, -23, -23, -23, -23, 
-23, -23, -23, -23, -23, -23, -22, -22, -22, -22, -22, -22, -22, 
-22, -22, -22, -22, -21, -21, -21, -21, -21, -21, -21, -21, -21, 
-21, -21, -20, -20, -20, -20, -20, -20, -20, -20, -20, -20, -20, 
-19, -19, -19, -19, -19, -19, -19, -19, -19, -19, -19, -18, -18, 
-18, -18, -18, -18, -18, -18, -18, -18, -18, -17, -17, -17, -17, 
-17, -17, -17, -17, -17, -17, -17, -16, -16, -16, -16, -16, -16, 
-16, -16, -16, -16, -16, -15, -15, -15, -15, -15, -15, -15, -15, 
-15, -15, -15, -14, -14, -14, -14, -14, -14, -14, -14, -14, -14, 
-14, -13, -13, -13, -13, -13, -13, -13, -13, -13, -13, -13, -12, 
-12, -12, -12, -12, -12, -12, -12, -12, -12, -12, -11, -11, -11, 
-11, -11, -11, -11, -11, -11, -11, -11, -10, -10, -10, -10, -10, 
-10, -10, -10, -10, -10, -10, -9, -9, -9, -9, -9, -9, -9, -9, 
-9, -9, -9, -8, -8, -8, -8, -8, -8, -8, -8, -8, -8, -8, -7, -7, 
-7, -7, -7, -7, -7, -7, -7, -7, -7, -6, -6, -6, -6, -6, -6, -6, 
-6, -6, -6, -6, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -4, 
-4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -3, -3, -3, -3, -3, -3, 
-3, -3, -3, -3, -3, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, 
-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0), confidence = c(0, 0.8, 0.916666666666667, 1, 
1, 0.8, NaN, 1, NaN, NaN, NaN, 0, 0.756756756756757, 0.9, 1, 
1, 0.857142857142857, NaN, 1, NaN, NaN, NaN, 0, 0.756756756756757, 
0.9, 1, 1, 0.857142857142857, NaN, 1, NaN, NaN, NaN, 0, 0.787878787878788, 
0.785714285714286, 1, 1, 0.857142857142857, 1, 1, NaN, NaN, NaN, 
0, 0.774193548387097, 0.785714285714286, 1, 1, 0.8, 1, 1, NaN, 
NaN, NaN, 0, 0.774193548387097, 0.785714285714286, 1, 1, 0.8, 
1, 1, NaN, NaN, NaN, 0, 0.758620689655172, 0.785714285714286, 
1, 1, 0.8, 1, 1, NaN, NaN, NaN, 0, 0.758620689655172, 0.75, 1, 
1, 0.8, 1, 1, 1, NaN, NaN, 0, 0.740740740740741, 0.6875, 1, 1, 
0.888888888888889, 1, 1, 1, NaN, NaN, 0, 0.75, 0.5625, 1, 1, 
0.8, 1, 1, NaN, 1, NaN, 0, 0.695652173913044, 0.642857142857143, 
0.875, 1, 0.8, 1, 1, NaN, 1, NaN, 0, 0.666666666666667, 0.736842105263158, 
0.846153846153846, 0.8, 0.8, 1, 1, NaN, 1, NaN, 0, 0.615384615384616, 
0.736842105263158, 0.8, 0.777777777777778, 0.909090909090909, 
1, 1, NaN, 1, 1, 0, 0.615384615384616, 0.736842105263158, 0.8, 
0.777777777777778, 0.8, 1, 1, 1, 1, 1, 0, 0.6, 0.727272727272727, 
1, 0.636363636363636, 0.8, 1, 1, 1, 1, 1, 0, 0.333333333333333, 
0.769230769230769, 1, 0.555555555555555, 0.857142857142857, 1, 
1, 1, 1, 1, 0, 0.5, 0.666666666666667, 0.823529411764706, 0.555555555555555, 
0.8, 1, 1, 1, 1, 1, 0, 0.5, 0.615384615384616, 1, 0.5625, 0.888888888888889, 
1, 1, 1, 1, 1, 0, 0.5, 0.444444444444445, 1, 0.615384615384616, 
0.785714285714286, 1, NaN, 1, 1, 1, 0, 0.5, 0.5, 0.857142857142857, 
0.666666666666667, 0.583333333333333, 1, 1, 1, 1, 1, 0, 0, 0.5, 
0.8, 0.727272727272727, 0.615384615384616, 1, 1, 1, 1, 1, 0, 
0, 0.5, 0.571428571428572, 0.769230769230769, 0.545454545454546, 
1, 1, 1, 1, NaN, 0, 0, 0, 0, 0.769230769230769, 0.666666666666667, 
0.5, 1, 1, 1, 1, 0, 0, 0, 0, 0.666666666666667, 0.666666666666667, 
NaN, 0.818181818181818, 1, 1, 1, 0, 0, 0, 0, 0, 0.666666666666667, 
1, 0.666666666666667, 1, 1, 1, 0, 0, 0, 0, 0, 0.363636363636364, 
1, 0.5, 1, 1, 1, 0, 0, 0, 0, 0, 0.222222222222222, 1, 0.5, 1, 
1, 1, 0, 0, NaN, 0, 0, 0.222222222222222, NaN, 0.5, 1, 1, 1, 
0, 0, NaN, 0, 0, 0, NaN, 0, 0.8, NaN, 1, 0, 0, NaN, 0, 0, 0, 
NaN, NaN, 0.5, 1, NaN, 0, 0, NaN, 0, 0, 0, NaN, NaN, 0, 0.666666666666667, 
NaN)), .Names = c("counts", "time", "confidence"), row.names = c(NA, 
-341L), class = "data.frame")

# before smoothing, the graph is really ugly

ct <- .2 # color threshold
before <- ggplot(df, aes(time, counts, fill=confidence)) + 
				geom_tile() +
				scale_fill_gradientn("Confidence", colours=c('darkblue', 'lightblue', 'darkred', 'red'), 
				    values=c(1, ct+.001,ct,0), rescale=FALSE, 
				    breaks=c(0, ct/2, ct, ct*1.5, ct*2, 1)) +
				scale_x_continuous("Days", breaks=(-4:0)*7) +
				scale_y_continuous("Counts") +
				opts(title='Before')

# two smoothing parameters:
# Smoothing Nearest Neighbor parameter (larger is smoother)
smooth.nn <- .3
# Smoothing monotone bandwidth
smooth.bw <- .3

# use locfit to generate a kernel-based smoothing

smoothed.fit <- locfit(confidence ~ lp(counts, time, nn=smooth.nn), data=df)
smoothed.confidence <- predict(smoothed.fit, subset(df, select=c('counts', 'time')))
df2 <- df
df2$confidence <- smoothed.confidence

smoothed <- before %+% df2 + opts(title='Smoothed')

# use monoproc to force the smoothed surface to be monotonic

smoothed.fit.mono <- monoproc(list(counts=0:10, time=-30:0, confidence=smoothed.confidence), 
                          bandwidth=smooth.bw, mono1='increasing', mono2='decreasing', dir='xy')
df3 <- df2
df3$confidence <- pmin(1, pmax(0, smoothed.fit.mono@fit@z))

mono.smoothed <- smoothed %+% df3 + opts(title='Monotonic + Smoothed')

cat('to display the graphs:
print(before)
print(smoothed)
print(mono.smoothed)')