[R] add median value and standard deviation bar to lattice plot

Bert Gunter bgunter.4567 at gmail.com
Fri Mar 17 15:59:08 CET 2017


If I understand you correcty, I suggest you consult a local
statistician. Individual data points *cannot* have "confidence
intervals."  So you clearly need some statistical guidance, which is
not what this list is about.

Cheers,
Bert


Bert Gunter

"The trouble with having an open mind is that people keep coming along
and sticking things into it."
-- Opus (aka Berkeley Breathed in his "Bloom County" comic strip )


On Fri, Mar 17, 2017 at 1:27 AM, Luigi Marongiu
<marongiu.luigi at gmail.com> wrote:
> dear Bert,
> the code as I read it draws 30 median lines, one for each cluster.
> what I am looking for is a function that will plot 190 individual
> confidence intervals, one for each measurement. would the code cover
> even that instance?
> regards
> luigi
>
> On Thu, Mar 16, 2017 at 2:41 PM, Bert Gunter <bgunter.4567 at gmail.com> wrote:
>> Just add whatever further code to decorate the groups as you like
>> within the panel.groups function. I believe I have given you
>> sufficient information in my code for you to do that if you study the
>> code carefully. Depending on what you decide to do -- which is
>> statistical and OT here (and not something I would offer specific
>> advice on remotely anyway) -- you may also have to pass down
>> additional arguments based on computations that you do with *all* the
>> data from *all* groups together.
>>
>> Cheers,
>> Bert
>>
>>
>> Bert Gunter
>>
>> "The trouble with having an open mind is that people keep coming along
>> and sticking things into it."
>> -- Opus (aka Berkeley Breathed in his "Bloom County" comic strip )
>>
>>
>> On Thu, Mar 16, 2017 at 1:38 AM, Luigi Marongiu
>> <marongiu.luigi at gmail.com> wrote:
>>> dear Bert,
>>> thank you for the solution, it worked perfectly. However I still would
>>> like to know how reliable are the dots that are plotted, that is why i
>>> would like to have individual bars on each dot (if possible). the
>>> standard deviation maybe is not the right tool and the confidence
>>> interval is perhaps better, but the procedure should be the same: draw
>>> an arrow from the lower to the upper limit. is that possible?
>>> regards,
>>> luigi
>>>
>>> PS sorry for the formatting, usually plain text is my default; it
>>> should have switched to html when i replied to a previous email but
>>> the difference does not show up when i type...
>>>
>>> On Wed, Mar 15, 2017 at 4:28 PM, Bert Gunter <bgunter.4567 at gmail.com> wrote:
>>>> There may be a specific function that handles this for you, but to
>>>> roll your own, you need a custom panel.groups function, not the
>>>> default. You need to modify the panel function (which is
>>>> panel.superpose by default) to pass down the "col" argument to the
>>>> panel.segments call in the panel.groups function.
>>>>
>>>> This should get you started:
>>>>
>>>> useOuterStrips(
>>>>    strip = strip.custom(par.strip.text = list(cex = 0.75)),
>>>>    strip.left = strip.custom(par.strip.text = list(cex = 0.75)),
>>>>    stripplot(
>>>>       average ~ type|target+cluster,
>>>>       panel = function(x,y,col,...)
>>>>          panel.superpose(x,y,col=col,...),
>>>>       panel.groups = function(x,y,col,...){
>>>>          panel.stripplot(x,y,col=col,...)
>>>>          m <- median(y)
>>>>          panel.segments(x0 = x[1] -.5, y0 = m,
>>>>                         x1 = x[1] +.5, y1 = m,
>>>>                         col=col, lwd=2
>>>>                         )
>>>>       },
>>>>       my.data,
>>>>       groups = type,
>>>>       pch=1,
>>>>       jitter.data = TRUE,
>>>>       main = "Group-wise",
>>>>       xlab = expression(bold("Target")), ylab = expression(bold("Reading")),
>>>>       col = c("grey", "green", "red"),
>>>>       par.settings = list(strip.background = list(col=c("paleturquoise",
>>>>                                                         "grey"))),
>>>>       scales = list(alternating = FALSE, x=list(draw=FALSE)),
>>>>       key = list(
>>>>          space = "top",
>>>>          columns = 3,
>>>>          text = list(c("Blank", "Negative", "Positive"), col="black"),
>>>>          rectangles = list(col=c("grey", "green", "red"))
>>>>       )
>>>>    )
>>>> )
>>>>
>>>> FWIW, I think adding 1 sd bars is a bad idea statistically.
>>>>
>>>> And though it made no difference here, please post in pain text, not HTML.
>>>>
>>>> Bert Gunter
>>>>
>>>> "The trouble with having an open mind is that people keep coming along
>>>> and sticking things into it."
>>>> -- Opus (aka Berkeley Breathed in his "Bloom County" comic strip )
>>>>
>>>>
>>>> On Wed, Mar 15, 2017 at 2:22 AM, Luigi Marongiu
>>>> <marongiu.luigi at gmail.com> wrote:
>>>>> Dear all,
>>>>> I am analyzing some multivariate data that is organized like this:
>>>>> 1st variable = cluster (A or B)
>>>>> 2nd variable = target (a, b, c, d, e)
>>>>> 3rd variable = type (blank, negative, positive)
>>>>> 4th variable = sample (the actual name of the sample)
>>>>> 5th variable = average (the actual reading -- please not that this is the
>>>>> mean of different measures with an assumed normal distribution, but the
>>>>> assumption might not always be true)
>>>>> 6th variable = stdev (the standard deviation associated with each reading)
>>>>> 7th variable = ll (lower limit that is average stdev)
>>>>> 8th variable = ul (upper limit that is average + stdev)
>>>>>
>>>>> I am plotting the data using lattice's stripplot and I would need to add:
>>>>> 1. an error bar for each measurement. the bar should be possibly coloured
>>>>> in light grey and semitransparent to reduce the noise of the plot.
>>>>> 2. a type-based median bar to show differences in measurements between
>>>>> blanks, negative and positive samples within each panel.
>>>>>
>>>>> How would I do that?
>>>>> Many thanks,
>>>>> Luigi
>>>>>
>>>>>>>>
>>>>> cluster <- c(rep("A", 90), rep("B", 100))
>>>>> sample <- c(
>>>>>   rep(c("cow-01", "cow-02", "cow-03", "cow-04", "cow-05", "cow-06",
>>>>> "cow-07", "cow-08", "cow-09", "cow-10", "cow-11",
>>>>>         "cow-12", "cow-13", "cow-14", "cow-15", "cow-16", "cow-17",
>>>>> "blank"), 5),
>>>>>   rep(c("cow-26", "cow-35", "cow-36", "cow-37", "cow-38", "cow-39",
>>>>> "cow-40", "cow-41", "cow-42", "cow-43", "cow-44", "cow-45",
>>>>>         "cow-46", "cow-47", "cow-48", "cow-49", "cow-50", "cow-51",
>>>>> "cow-59", "blank"), 5)
>>>>> )
>>>>> type <- c(
>>>>>   rep(c("negative", "negative", "negative", "negative", "negative",
>>>>> "negative", "negative", "negative", "positive", "positive",
>>>>>         "positive", "positive", "positive", "positive", "positive",
>>>>> "positive", "positive", "blank"), 5),
>>>>>   rep(c("negative", "positive", "negative", "negative", "negative",
>>>>> "negative", "negative", "negative", "positive", "positive",
>>>>>         "positive", "positive", "positive", "positive", "positive",
>>>>> "positive", "positive", "positive", "positive", "blank"), 5)
>>>>> )
>>>>> target <- c(
>>>>> c(rep("a", 18), rep("b", 18), rep("c", 18), rep("d", 18), rep("e", 18)),
>>>>> c(rep("a", 20), rep("b", 20), rep("c", 20), rep("d", 20), rep("e", 20))
>>>>> )
>>>>> average <- c(88.5, 49, 41, 33, 35, 45, 95, 30, 41, 64, 22, 29, 59, 71, 128, 39,
>>>>> 42, 47, 86, 100,
>>>>>              69, 44, 53, 66, 66, 71, 161, 69, 22.5, 30, 67, 99, 129, 94, 49,
>>>>> 33, 28, 31, 26, 23,
>>>>>              30, 41, 35, 23, 38, 43, 15, 21, 45, 51.5, 34, 26, 43, 32.5, 59,
>>>>> 58.5, 61, 62.5, 58,
>>>>>              59.5, 60.5, 60, 64, 110, 55, 66, 197, 83.5, 155, 76, 125, 90, 73,
>>>>> 84, 95.5, 62, 82, 138,
>>>>>              103.5, 57, 138, 149.5, 57, 54, 245.5, 191, 131, 96, 176, 45, 76,
>>>>> 33, 37, 51, 44, 50, 54,
>>>>>              66, 49, 90, 66.5, 42.5, 67, 56, 54, 50, 45, 99, 50, 51.5, 212, 40,
>>>>> 68, 121, 80, 57,
>>>>>              81.5, 128, 77, 119.5, 126, 184, 101, 103, 88, 100, 140, 186, 297,
>>>>> 32, 184, 36, 45, 45, 44,
>>>>>              86, 65, 61, 76, 62, 136, 84, 80, 56, 109, 116, 54, 59,
>>>>> 79, 34, 74.5,
>>>>> 54, 49, 55, 56,
>>>>>              59, 56, 56, 57, 67, 65, 63, 52, 58, 59, 56, 54, 66, 92, 87, 59,
>>>>> 33, 58, 51, 54,
>>>>>              52, 47, 45, 42, 52, 57, 79, 42, 45.5, 47, 47, 36, 50, 53, 49 )
>>>>> stdev <- c(17.85, 6.31, 3.42, 1.04, 0.51, 6.04, 38.43, 2.78, 5.55, 26.72, 1.83,
>>>>> 9.92, 4.59, 19, 7.96,
>>>>>                7.5, 1.06, 9.66, 75.94, 36.79, 50.45, 9.79, 1.55, 11.42, 64.12,
>>>>> 0.79, 15.14, 16.15, 8.12, 4.04, 92.57, 35.35,
>>>>>                42.28, 52.96, 7.06, 4.97, 1.15, 4.77, 6.59, 7.27, 0.75, 4.25,
>>>>> 9, 0.1, 1.14, 4.17, 6.73, 3.81, 3.27,
>>>>>                97.44, 9.74, 0.45, 8.14, 5.91, 13.1, 98.22, 8.92, 72.62, 70.26,
>>>>> 59.46, 29.89, 56.35, 91.25, 49.94, 20.65, 62.04,
>>>>>                95.13, 35.89, 99.64, 29.44, 33.12, 45.91, 96.69, 9.05, 38.56,
>>>>> 3.09, 0.6, 8.69, 16.95, 74.03, 84.05, 39.87, 15.52,
>>>>>                27.92, 35.72, 80.26, 71.93, 66.73, 87.8, 5.43, 98.3, 7.41, 9.86,
>>>>> 63.64, 0.36, 5.84, 1.58, 20.1, 4.21, 82.12,
>>>>>                19.29, 9.02, 22.12, 54.08, 74.95, 3.24, 9.67, 67.98,
>>>>> 9.92, 40.69,
>>>>> 6.24, 8.76, 74.25, 46.34, 25.69, 90.63, 83.71,
>>>>>                73.53, 57.88, 15.84, 82.07, 67.45, 47.39, 98.77, 75.1,
>>>>> 64.9, 3.71,
>>>>> 87.44, 61.06, 4.77, 57.54, 7.68, 4.54, 6.15,
>>>>>                3.32, 60.39, 33.78, 66.22, 18.67, 76.53, 63.54, 47.06, 38.47,
>>>>> 88.15, 18.25, 4.26, 67.19, 88.87, 29.65, 7.33, 68.18,
>>>>>                28.03, 6.91, 77.82, 22.23, 73.23, 95.21, 27.11, 37.01, 34.88,
>>>>> 28.15, 11.27, 15.67, 96.08, 89.52, 28.6, 8.22, 23.55,
>>>>>                59.2, 36.38, 41.38, 0.4, 56.82, 32.35, 20.6, 18.13, 8.15, 1.08,
>>>>> 9.85, 1.07, 37.75, 97.6, 7.16, 8.51, 4.42,
>>>>>                0.15, 1.28, 7.42, 71.15, 9.39)
>>>>> ll <- c(70.65, 42.69, 37.58, 31.96, 34.49, 38.96, 56.57, 27.22, 35.45,
>>>>> 37.28, 20.17, 19.08, 54.41, 52, 120.04, 31.5, 40.94, 37.34,
>>>>>         10.06, 63.21, 18.55, 34.21, 51.45, 54.58, 1.88, 70.21, 145.86,
>>>>> 52.85, 14.38, 25.96, -25.57, 63.65, 86.72, 41.04, 41.94, 28.03,
>>>>>         26.85, 26.23, 19.41, 15.73, 29.25, 36.75, 26, 22.9, 36.86, 38.83,
>>>>> 8.27, 17.19, 41.73, -45.94, 24.26, 25.55, 34.86, 26.59, 45.9,
>>>>>         -39.72, 52.08, -10.12, -12.26, 0.0399999999999991, 30.61, 3.65,
>>>>> -27.25, 60.06, 34.35, 3.96, 101.87, 47.61, 55.36, 46.56, 91.88, 44.09,
>>>>>         -23.69, 74.95, 56.94, 58.91, 81.4, 129.31, 86.55, -17.03, 53.95,
>>>>> 109.63, 41.48, 26.08, 209.78, 110.74, 59.07, 29.27, 88.2, 39.57,
>>>>>         -22.3, 25.59, 27.14, -12.64, 43.64, 44.16, 52.42, 45.9, 44.79, 7.88,
>>>>> 47.21, 33.48, 44.88, 1.92, -20.95, 46.76, 35.33, 31.02,
>>>>>         40.08, 10.81, 205.76, 31.24, -6.25, 74.66, 54.31, -33.63,
>>>>> -2.20999999999999, 54.47, 19.12, 103.66, 43.93, 116.55, 53.61, 4.23,
>>>>>         12.9, 35.1, 136.29, 98.56, 235.94, 27.23, 126.46, 28.32, 40.46,
>>>>> 38.85, 40.68, 25.61, 31.22, -5.22, 57.33, -14.53, 72.46, 36.94,
>>>>>         41.53, -32.15, 90.75, 111.74, -13.19, -29.87, 49.35, 26.67,
>>>>> 6.31999999999999, 25.97, 42.09, -22.82, 33.77, -14.23, -39.21, 28.89,
>>>>>         19.99, 32.12, 36.85, 51.73, 36.33, -38.08, -30.52, 27.4, 45.78,
>>>>> 42.45, 32.8, 50.62, 17.62, 32.6, 1.18, 18.65, 33.4, 33.87, 38.85,
>>>>>         43.92, 32.15, 50.93, 19.25, -18.6, 34.84, 36.99, 42.58, 46.85,
>>>>> 34.72, 42.58, -18.15, 39.61)
>>>>> ul <- c(106.35, 55.31, 44.42, 34.04, 35.51, 51.04, 133.43, 32.78, 46.55,
>>>>> 90.72, 23.83, 38.92, 63.59, 90, 135.96, 46.5, 43.06, 56.66,
>>>>>         161.94, 136.79, 119.45, 53.79, 54.55, 77.42, 130.12, 71.79, 176.14,
>>>>> 85.15, 30.62, 34.04, 159.57, 134.35, 171.28, 146.96, 56.06, 37.97,
>>>>>         29.15, 35.77, 32.59, 30.27, 30.75, 45.25, 44, 23.1, 39.14, 47.17,
>>>>> 21.73, 24.81, 48.27, 148.94, 43.74, 26.45, 51.14, 38.41, 72.1,
>>>>>         156.72, 69.92, 135.12, 128.26, 118.96, 90.39, 116.35, 155.25,
>>>>> 159.94, 75.65, 128.04, 292.13, 119.39, 254.64, 105.44, 158.12, 135.91,
>>>>> 169.69,
>>>>>         93.05, 134.06, 65.09, 82.6, 146.69, 120.45, 131.03, 222.05, 189.37,
>>>>> 72.52, 81.92, 281.22, 271.26, 202.93, 162.73, 263.8, 50.43, 174.3,
>>>>>         40.41, 46.86, 114.64, 44.36, 55.84, 55.58, 86.1, 53.21, 172.12,
>>>>> 85.79, 51.52, 89.12, 110.08, 128.95, 53.24, 54.67, 166.98, 59.92,
>>>>>         92.19, 218.24, 48.76, 142.25, 167.34, 105.69, 147.63, 165.21,
>>>>> 201.53, 134.88, 135.34, 208.07, 251.45, 148.39, 201.77, 163.1, 164.9,
>>>>> 143.71,
>>>>>         273.44, 358.06, 36.77, 241.54, 43.68, 49.54, 51.15, 47.32, 146.39,
>>>>> 98.78, 127.22, 94.67, 138.53, 199.54, 131.06, 118.47, 144.15, 127.25,
>>>>>         120.26, 121.19, 147.87, 108.65, 41.33, 142.68, 82.03, 55.91, 132.82,
>>>>> 78.23, 132.23, 151.21, 83.11, 94.01, 101.88, 93.15, 74.27, 67.67,
>>>>>         154.08, 148.52, 84.6, 62.22, 89.55, 151.2, 123.38, 100.38, 33.4,
>>>>> 114.82, 83.35, 74.6, 70.13, 55.15, 46.08, 51.85, 53.07, 94.75, 176.6,
>>>>>         49.16, 54.01, 51.42, 47.15, 37.28, 57.42, 124.15, 58.39)
>>>>> my.data <- data.frame(cluster, type, target, sample, average, stdev, ll,
>>>>> ul, stringsAsFactors = FALSE)
>>>>>
>>>>> library(lattice)
>>>>> library(latticeExtra)
>>>>> useOuterStrips(
>>>>>   strip = strip.custom(par.strip.text = list(cex = 0.75)),
>>>>>   strip.left = strip.custom(par.strip.text = list(cex = 0.75)),
>>>>>   stripplot(
>>>>>     average ~ type|target+cluster,
>>>>>     my.data,
>>>>>     groups = type,
>>>>>     pch=1,
>>>>>     jitter.data = TRUE,
>>>>>     main = "Group-wise",
>>>>>     xlab = expression(bold("Target")), ylab = expression(bold("Reading")),
>>>>>     col = c("grey", "green", "red"),
>>>>>     par.settings = list(strip.background = list(col=c("paleturquoise",
>>>>> "grey"))),
>>>>>     scales = list(alternating = FALSE, x=list(draw=FALSE)),
>>>>>     key = list(
>>>>>       space = "top",
>>>>>       columns = 3,
>>>>>       text = list(c("Blank", "Negative", "Positive"), col="black"),
>>>>>       rectangles = list(col=c("grey", "green", "red"))
>>>>>     )
>>>>>   )
>>>>> )
>>>>>
>>>>>         [[alternative HTML version deleted]]
>>>>>
>>>>> ______________________________________________
>>>>> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
>>>>> https://stat.ethz.ch/mailman/listinfo/r-help
>>>>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>>>>> and provide commented, minimal, self-contained, reproducible code.



More information about the R-help mailing list