# Comparison method for curves using percentiles

I have a question regarding the use of quantiles for determining the enveloppe of a curve. This is what I am doing: I have a continuous variable ("var") that I turned discrete using cut ("var_cut"), and a related variable obtained from the the continuous variable ("modvar"). What I'm doing is plotting modvar~var_cut, and I want to have a sense of the variability of modvar. Here is the method I chose:

#df containing var and modvar
structure(list(var = c(0.1968, 0.2263667, 0.1769, 0.2318, 0.2001333,
0.2382667, 0.2005, 0.2022667, 0.1699333, 0.2115667, 0.212, 0.2218667,
0.2327333, 0.2224333, 0.1690333, 0.1961333, 0.1756667, 0.2268333,
0.1938667, 0.1983, 0.1914333, 0.1745333, 0.2382, 0.2068333, 0.2509333,
0.221, 0.2075667, 0.2475333, 0.2463333, 0.2354, 0.2335, 0.2382,
0.2636667, 0.1829667, 0.2180333, 0.1703333, 0.2177333, 0.1932667,
0.2281, 0.1960667, 0.1975333, 0.1640333, 0.2021667, 0.2044333,
0.2124, 0.2267, 0.2202333, 0.1648667, 0.1898, 0.168, 0.2225,
0.1899667, 0.1966667, 0.183, 0.1678667, 0.2288333, 0.2006, 0.2389333,
0.2105, 0.2018667, 0.2457667, 0.2393333, 0.2286, 0.2280333, 0.2319,
0.2565333, 0.1838, 0.2189667, 0.1710667, 0.2184, 0.194, 0.2289333,
0.1968, 0.1984, 0.1646667, 0.2029667, 0.2053667, 0.2132333, 0.2274667,
0.2211, 0.1655333, 0.1907333, 0.1688333, 0.2234, 0.1908, 0.1975333,
0.1838333, 0.1686, 0.2297333, 0.2013667, 0.2397667, 0.2113333,
0.2027333, 0.2467333, 0.2402, 0.2295333, 0.2289333, 0.2328333,
0.2574333, 0.1795667), modvar = c(1.01575728698598, 0.978902741156023,
1.04056240429755, 0.972130196236979, 1.01160236751187, 0.964069530301364,
1.01114528024965, 1.00894310935747, 1.04924631438672, 0.997350768101313,
0.99681066471784, 0.984511938538037, 0.97096684869995, 0.983805678263226,
1.05036815386312, 1.01658832074033, 1.04209969832671, 0.978321129711924,
1.01941361113723, 1.0138875545253, 1.02244681578377, 1.04351246817399,
0.964152671071449, 1.00325089585421, 0.948280761510472, 0.985592269953813,
1.00233672132977, 0.95251882175466, 0.954014607723197, 0.967642838331369,
0.970011166114885, 0.964152671071449, 0.932408727300664, 1.03300033368478,
0.989290226814528, 1.04874771906387, 0.989664173306662, 1.0201615041215,
0.976742202973302, 1.01667133686158, 1.01484323711037, 1.05660059539869,
1.00906775818819, 1.00624246779128, 0.996312069394995, 0.978487286603262,
0.986547952538877, 1.05556177204354, 1.02448270513577, 1.05165615023086,
0.983722537493141, 1.02427491553498, 1.01592344387731, 1.03295882562415,
1.0518223071222, 0.975828153097695, 1.01102063141894, 0.963238621195842,
0.998680397178511, 1.00944170468032, 0.954720867998008, 0.962740025872996,
0.976118958819745, 0.976825343743386, 0.972005547406268, 0.941300426990633,
1.03196163497846, 0.988126754628668, 1.04783354453944, 0.988833139552309,
1.0192474542459, 0.975703504266984, 1.01575728698598, 1.01376290569459,
1.05581106970497, 1.00807056754249, 1.00507899560542, 0.995273370688676,
0.977531604018197, 0.985467621123101, 1.05473086293802, 1.02331935759875,
1.05061745152455, 0.982600698016739, 1.02323621682866, 1.01484323711037,
1.03192012691783, 1.0509082572466, 0.974706313621292, 1.01006494883388,
0.962199797840693, 0.997641698472193, 1.00836149791338, 0.953516012400351,
0.96165969445722, 0.974955611282715, 0.975703504266984, 0.970842199869238,
0.94017858751423, 1.03723839392897)), .Names = c("var", "modvar"
), row.names = c(NA, 100L), class = "data.frame")

#Calculation of discrete variable, as well as lower and upper boundaries of modvar
df$var_cut<-cut(df$var, quantile(df$var, (0:10)/10), include.lowest=TRUE) df$var_cut<-cut(df$var, quantile(df$var, (0:10)/10), include.lowest=TRUE, labels=c(1:length(levels(df$var_cut)))) df$lowervar<-ifelse(df$var_cut=="1",df$lowervar<-quantile(df[df$var_cut=="1","modvar"],c(0.05), na.rm=T), ifelse(df$var_cut=="2",df$lowervar<-quantile(df[df$var_cut=="2","modvar"],c(0.05), na.rm=T),
ifelse(df$var_cut=="3",df$lowervar<-quantile(df[df$var_cut=="3","modvar"],c(0.05), na.rm=T), ifelse(df$var_cut=="4",df$lowervar<-quantile(df[df$var_cut=="4","modvar"],c(0.05), na.rm=T),
ifelse(df$var_cut=="5",df$lowervar<-quantile(df[df$var_cut=="5","modvar"],c(0.05), na.rm=T), ifelse(df$var_cut=="6",df$lowervar<-quantile(df[df$var_cut=="6","modvar"],c(0.05), na.rm=T),
ifelse(df$var_cut=="7",df$lowervar<-quantile(df[df$var_cut=="7","modvar"],c(0.05), na.rm=T), ifelse(df$var_cut=="8",df$lowervar<-quantile(df[df$var_cut=="8","modvar"],c(0.05), na.rm=T),
ifelse(df$var_cut=="9",df$lowervar<-quantile(df[df$var_cut=="9","modvar"],c(0.05), na.rm=T), ifelse(df$var_cut=="10",df$lowervar<-quantile(df[df$var_cut=="10","modvar"],c(0.05), na.rm=T),NA))))))))))

df$uppervar<-ifelse(df$var_cut=="1",df$uppervar<-quantile(df[df$var_cut=="1","modvar"],c(0.95), na.rm=T),
ifelse(df$var_cut=="2",df$uppervar<-quantile(df[df$var_cut=="2","modvar"],c(0.95), na.rm=T), ifelse(df$var_cut=="3",df$uppervar<-quantile(df[df$var_cut=="3","modvar"],c(0.95), na.rm=T),
ifelse(df$var_cut=="4",df$uppervar<-quantile(df[df$var_cut=="4","modvar"],c(0.95), na.rm=T), ifelse(df$var_cut=="5",df$uppervar<-quantile(df[df$var_cut=="5","modvar"],c(0.95), na.rm=T),
ifelse(df$var_cut=="6",df$uppervar<-quantile(df[df$var_cut=="6","modvar"],c(0.95), na.rm=T), ifelse(df$var_cut=="7",df$uppervar<-quantile(df[df$var_cut=="7","modvar"],c(0.95), na.rm=T),
ifelse(df$var_cut=="8",df$uppervar<-quantile(df[df$var_cut=="8","modvar"],c(0.95), na.rm=T), ifelse(df$var_cut=="9",df$uppervar<-quantile(df[df$var_cut=="9","modvar"],c(0.95), na.rm=T),
ifelse(df$var_cut=="10",df$uppervar<-quantile(df[df$var_cut=="10","modvar"],c(0.95), na.rm=T),NA))))))))))  I have been working under the assumption that by using the "probs" argument of the quantile function, I can obtain an enveloppe for my modvar curve that is similar to a confidence interval, specifying the lower and upper bounds to correspond respectively to 0.05 and 0.95 probabilities in the quantile function. Would you say this is acceptable a method? What about comparing two different curves this way? I would like to check for overlaps between different variables for example. What I would be doing is plotting the same modvar~var_cut but for different databases, and checking if the lower and upper bounds of the respective curves overlap. Edit: adding the plot: png();with( df, matplot(x=var_cut, y=df[, c(lowervar, uppervar)] ) ); dev.off()  Edit: data and code for the graphs Here is the data for the other curve: structure(list(var = c(0.1949667, 0.2248, 0.1982, 0.2285, 0.2143, 0.2211, 0.2169333, 0.2341, 0.1923333, 0.2236333, 0.2548, 0.2683333, 0.2956, 0.2730667, 0.2101, 0.2247333, 0.1999, 0.2598, 0.2282333, 0.2730333, 0.2576, 0.2164667, 0.2526333, 0.2494667, 0.2534, 0.2212333, 0.2391333, 0.2781, 0.2585667, 0.1924667, 0.221, 0.1933, 0.2213, 0.2128, 0.2225333, 0.2110333, 0.2229, 0.1858667, 0.2215333, 0.2561, 0.257, 0.2855, 0.2605, 0.1991667, 0.2201333, 0.1917667, 0.2508, 0.2224667, 0.2560333, 0.2445667, 0.2066333, 0.2493333, 0.2407, 0.2503333, 0.2195, 0.2225667, 0.2782333, 0.2544667, 0.2097333, 0.2266, 0.1834333, 0.2147333, 0.2058667, 0.2255333, 0.2095, 0.2103, 0.1803667, 0.2125667, 0.2397, 0.2451333, 0.2535333, 0.2355333, 0.1825667, 0.2147667, 0.1873667, 0.2442667, 0.2135333, 0.2208333, 0.2074, 0.1857667, 0.2467, 0.2207667, 0.2559333, 0.2185, 0.2176333, 0.2641333, 0.2520333, 0.2162, 0.2329333, 0.1889667, 0.2209667, 0.2119667, 0.2316, 0.2155333, 0.2157, 0.1849667, 0.2184, 0.2456333, 0.2507667, 0.2601), modvar = c(1.04309748897633, 1.00608866182958, 1.03908651329433, 1.00149873512207, 1.01911412951303, 1.01067858853708, 1.0158474662701, 0.99455181902423, 1.04636427627134, 1.00753597736732, 0.968873039876843, 0.952084700657745, 0.918259794021125, 0.946212819826043, 1.02432431658642, 1.00617140456239, 1.03697762805035, 0.962670436218054, 1.00182958200123, 0.946254253218484, 0.965399581827921, 1.01642629324353, 0.971560876146343, 0.975489109095527, 0.970609768901304, 1.01051322712354, 0.988307906025073, 0.939968906826886, 0.964200370436531, 1.04619879080572, 1.01080264061026, 1.04516506487995, 1.01043048439073, 1.02097491061067, 1.00890055017225, 1.02316653858747, 1.00844565121992, 1.05438622763532, 1.01014107090401, 0.967260362925558, 0.966143894266976, 0.930789053411879, 0.961802071705824, 1.03788730190294, 1.01187779992847, 1.04706715531795, 0.973835122803874, 1.00898316885299, 0.967343105658366, 0.98156766068114, 1.0286248298072, 0.975654594561143, 0.986364382194628, 0.974414073829386, 1.01266342170789, 1.00885911677981, 0.939803545413343, 0.969286505436738, 1.02477921553875, 1.00385572451241, 1.05740491078398, 1.01857661187996, 1.02957581300017, 1.00517898797698, 1.02506862902547, 1.02407621244007, 1.06120909165999, 1.02126432409739, 0.987604902926386, 0.980864781634526, 0.970444407487761, 0.992773780659401, 1.05847994605012, 1.01853517848752, 1.05252544653769, 0.981939816900667, 1.02006523675807, 1.01100943541624, 1.02767372256216, 1.0545102797085, 0.978921257804081, 1.01109205409698, 0.967467157731542, 1.01390394243965, 1.01497910175787, 0.957294887731128, 0.972305188585397, 1.01675714012269, 0.995999134561971, 1.05054061336687, 1.01084394995062, 1.02200863653644, 0.997653120853624, 1.01758419529456, 1.01737740048857, 1.05550269629391, 1.01402799451283, 0.980244521268647, 0.973876432144242, 0.962298279998527)), .Names = c("var", "modvar"), row.names = c(NA, -100L), class = "data.frame")  So you can use the same code to obtain var_cut and the upper and lower boundaries. Here is the code I used for the graphs, modified so it could take into account the two different df: p<-ggplot(df1, aes(var_cut,modvar, group=1))+ geom_smooth(aes(color="red"), se=F, linetype="dotted", size=1)+ geom_line(data=df1,aes(var_cut,lowervar, color="red4"), size=1)+ geom_line(data=df1,aes(var_cut,uppervar, color="red4"), size=1)+ geom_ribbon(data=df1, aes(var_cut,ymin=lowervar,ymax=uppervar), fill="lightpink", alpha=0.4)+ geom_smooth(data=df2, aes(var_cut,modvar, group=1, color="green"), se=F, linetype="dotted", size=1)+ geom_line(data=df2,aes(var_cut,uppervar, color="green4"), size=1)+ geom_line(data=df2,aes(var_cut,lowervar, color="green4"), size=1)+ geom_ribbon(data=df2, aes(var_cut,ymin=lowervar,ymax=uppervar), fill="chartreuse1", alpha=0.4)+ ylim(min(df1$lowervar,df2$lowervar),max(df1$uppervar,df2\$uppervar))+
scale_colour_manual(name = 'Legend',
values =c('red'='red','green'='green', 'green4'='green4', 'red4'='red4'), labels = c('North','interval-N','South','interval-S'))+
scale_size_area() +
xlab("var by decile") +
ylab(expression(f[VPD.fall]))+
labs(title='modvar - North & South')


This is the code I would have used to do what I thought you were describing:

> aggregate(modvar ~ var_cut, df, function(x) quantile(x, c(0.05, 0.95)))
var_cut modvar.5% modvar.95%
1        1 1.0497511  1.0562453
2        2 1.0319388  1.0483363
3        3 1.0166257  1.0243892
4        4 1.0110767  1.0158487
5        5 1.0027481  1.0097845
6        6 0.9884446  0.9982130
7        7 0.9783959  0.9861179
8        8 0.9732894  0.9772138
9        9 0.9629644  0.9715381
10      10 0.9359052  0.9619568


Those are the 5th and 95%iles within each category of var_cut. That's a non-parameteric analogue of a 90% confidence interval. The plot looks like this:

with( grp.extremes, matplot(x=var_cut, y=grp.extremes[, c('modvar')] ) )


When interpreting such a result it is important to remember the sample sizes for these estimates. You only have 10 values per group, so estimating the 5th and 95th percentiles really is going way beyond what the data supports. There really is no datapoints on the other side of these values. You should probably have at least three or four values above and below those extremal thresholds for any stability in the estimates. With sample sizes of 10 per group I wouldn't think you could expect validity for much more than the 30th and 70th percentiles if there were any serious departure from "normality".

Your plot groups all of the quantile estimates for modvar at their var_cut value. The fact that they "go up" is kind of interesting, whereas you can see a fairly exact negative correlation when plotting:

 with(df, plot(var, modvar) )
> with(df, cor(var, modvar) )
[1] -1


You talk about "comparing two curves" but so far I do not have a good visualization of what these two curves really are. I've seen flawed regression analyses wehn reviewing papers for publication that reversed the apparent association of a relationship because the analyst created a grouping vector which was used as a predictor for the "outcome".

• Thanks for your reply. This is just a sample of the data I have at hand, and I have between 100 and 800 values per group depending on the database I am using, which I assume would be enough to use the 5th and 95th percentiles, am I wrong? Also, I'm interested in knowing why the values for the 5th and 95th percentiles are different using your method with aggregate. And is there any supporting litterature to legitimate the use of what you called "a non-parameteric analogue of a 90% confidence interval"? – Chris. Z Apr 3 '15 at 18:55
• I'll add some commentary on what I think you should have done to plot your results. The supporting literature would be found in texts on nonparametric statistics which would allow you to put confidence intervals around various quantiles and order statistics. A 90% CI in normal theory applications should give you the 5th and 95th percentiles if the data were perfectly normal. Any basic stats book should give you that. – DWin Apr 3 '15 at 19:14
• I will look into it; in the mean time, here is how I plotted a comparison for two databases, as I explained it in my post. This is fairly straight-forward and I don't know how I could improve on this but I'm open to your commentary! i.imgur.com/be88pQ8.png – Chris. Z Apr 3 '15 at 19:25
• Nicer looking plot that what I came up with, although It seems to be disguising the fact that these two variables are exactly correlated. – DWin Apr 3 '15 at 19:38
• If you had posted the code to create the plot it might have prevented confusion on my part. – DWin Apr 3 '15 at 19:51