Блог

Мы проводим урок статистики для студентов-биологов и пытаемся использовать R в качестве платформы для вычислений и визуализации данных. Насколько это возможно, мы хотели бы избегать использования дополнительных пакетов и делать что-либо ужасно «причудливое» в R; основное внимание в курсе уделяется статистике, а не программированию. Тем не менее, мы не нашли очень хорошего способа создания графика полосы ошибок в R для двухфакторного дизайна ANOVA. Для построения графика мы используем пакет ggplot2, и хотя в нем есть встроенный метод stat_summary для генерации 95% CI errorbar, способ их вычисления не всегда может быть правильным . Ниже я вручную просматриваю код для ANOVA и также вручную вычисляю 95% CI (со стандартной ошибкой, оцениваемой по общей остаточной дисперсии, а не только по внутригрупповой дисперсии, которую будет использовать сводный метод ggplot). В конце на самом деле есть график.

Итак, вопрос в том … есть ли более простой / быстрый / простой способ сделать все это?

 #   LIZARD LENGTH DATA
island.1 <- c(0.2, 5.9, 6.1, 6.5)
island.2 <- c(5.6, 14.8, 15.5, 16.4)
island.3 <- c(0.8, 3.9, 4.3, 4.9)
sex.codes <- c("Male", "Female", "Male", "Female")

#   PUTTING DATA TOGETHER IN A DATA FRAME
df.1 <- data.frame(island.1, island.2, island.3, sex.codes)

#   MELTING THE DATA FRAME INTO LONG FORM
library(reshape)
df.2 <- melt(df.1)

#   MEAN BY CELL
mean.island1.male <- with(df.2, mean(value[variable == "island.1" amp; sex.codes == "Male"]))
mean.island1.female <- with(df.2, mean(value[variable == "island.1" amp; sex.codes == "Female"]))
mean.island2.male <- with(df.2, mean(value[variable == "island.2" amp; sex.codes == "Male"]))
mean.island2.female <- with(df.2, mean(value[variable == "island.2" amp; sex.codes == "Female"]))
mean.island3.male <- with(df.2, mean(value[variable == "island.3" amp; sex.codes == "Male"]))
mean.island3.female <- with(df.2, mean(value[variable == "island.3" amp; sex.codes == "Female"]))

#   ADDING CELL MEANS TO DATA FRAME
df.2$means[df.2$variable == "island.1" amp; df.2$sex.codes == "Male"] <- mean.island1.male
df.2$means[df.2$variable == "island.1" amp; df.2$sex.codes == "Female"] <- mean.island1.female
df.2$means[df.2$variable == "island.2" amp; df.2$sex.codes == "Male"] <- mean.island2.male
df.2$means[df.2$variable == "island.2" amp; df.2$sex.codes == "Female"] <- mean.island2.female
df.2$means[df.2$variable == "island.3" amp; df.2$sex.codes == "Male"] <- mean.island3.male
df.2$means[df.2$variable == "island.3" amp; df.2$sex.codes == "Female"] <- mean.island3.female

#   LINEAR MODEL
lizard.model <- lm(value ~ variable*sex.codes, data=df.2)

#   CALCULATING RESIDUALS BY HAND:
df.2$residuals.1 <- df.2$value - df.2$means

#   CONFIRMING RESIDUALS FROM LINEAR MODEL:
df.2$residuals.2 <- residuals(lizard.model)

#   TWO FACTOR MAIN EFFECT ANOVA
lizard.anova <- anova(lizard.model)        

#   INTERACTION PLOT
interaction.plot(df.2$variable, df.2$sex.codes, df.2$value)

#   SAMPLE SIZE IN EACH CELL
n <- length(df.2$value[df.2$variable == "island.1" amp; df.2$sex.codes == "Male"])
# > n
# [1] 2

#   NOTE: JUST FOR CLARITY, PRETEND n=10
n <- 10

#   CALCULATING STANDARD ERROR
island.se <- sqrt(lizard.anova$M[4]/n)

#   HALF CONFIDENCE INTERVAL
island.ci.half <- qt(0.95, lizard.anova$D[4]) * island.se

#   MAKING SUMMARY DATA FRAME
summary.df <- data.frame(
        Means = c(mean.island1.male,
                mean.island1.female,
                mean.island2.male,
                mean.island2.female,
                mean.island3.male,
                mean.island3.female),
        Location = c("island1",
                "island1",
                "island2",
                "island2",
                "island3",
                "island3"),
        Sex = c("male",
                "female",
                "male",
                "female",
                "male",
                "female"),      
        CI.half = rep(island.ci.half, 6)        
        )

# > summary.df
# Means Location    Sex  CI.half
# 1  3.15  island1   male 2.165215
# 2  6.20  island1 female 2.165215
# 3 10.55  island2   male 2.165215
# 4 15.60  island2 female 2.165215
# 5  2.55  island3   male 2.165215
# 6  4.40  island3 female 2.165215

#   GENERATING THE ERRORBAR PLOT
library(ggplot2)

qplot(data=summary.df,
        y=Means,
        x=Location,
        group=Sex,
        ymin=Means-CI.half,
        ymax=Means CI.half,
        geom=c("point", "errorbar", "line"),
        color=Sex,
        shape=Sex,
        width=0.25)   theme_bw()
  

Вот еще одна попытка с использованием пакета sciplot. Альтернативные способы вычисления доверительных интервалов могут быть переданы в параметре ci.fun.

 lineplot.CI(variable,value, group =sex.codes , data = df.2, cex = 1.5,
            xlab = "Location", ylab = "means", cex.lab = 1.2, x.leg = 1,
            col = c("blue","red"), pch = c(16,16))
  

Я должен признать, что я совершенно сбит с толку вашим кодом. Не воспринимайте это как личную критику, но я настоятельно советую вам научить своих учеников максимально использовать возможности R. Они могут только выиграть от этого, и мой опыт показывает, что они быстрее понимают, что происходит, если я не забрасываю их головы строками и строками кода.

Прежде всего, вам не нужно вычислять средние вручную. Просто сделайте:

 df.2$mean <- with(df.2,ave(value,sex.codes,variable,FUN=mean))
  

Смотрите также ?ave . Это более понятно, чем беспорядок кода в вашем примере. Если у вас есть lizard.model, вы можете просто использовать

 fitted(lizard.model)
  

и сравните эти значения со средними.

Тогда я с вами категорически не согласен. То, что вы вычисляете, не является стандартной ошибкой вашего прогноза. Чтобы сделать это правильно, используйте predict() функцию

 outcome <- predict(lizard.model,se.fit=TRUE)
df.2$CI.half <- outcome$se / 2
  

Чтобы получить доверительный интервал для прогнозируемых средних, вы должны использовать правильные формулы, если хотите, чтобы ваши ученики правильно это поняли. Взгляните на раздел 3.5 этой невероятно отличной практической регрессии и Anova с использованием R из Faraway. Он содержит тонны примеров кода, где все вычисляется вручную удобным и кратким способом. Он будет полезен как вам, так и вашим ученикам. Я многому научился из этого и часто использую его в качестве руководства при объяснении этих вещей студентам.

Теперь, чтобы получить сводный фрейм данных, у вас есть несколько вариантов, но этот работает и вполне понятен.

 summary.df <- unique(df.2[,-c(3,5,6)])
names(summary.df) <- c('Sex','Location','Means','CI.half')
  

И теперь вы можете просто запустить свой код графика в том виде, в каком он есть.

В качестве альтернативы, если вам нужна ошибка прогнозирования для ваших значений, вы можете использовать следующее:

 lizard.predict <- predict(lizard.model,interval='prediction')
df.2$lower <- lizard.predict[,2]
df.2$upper <- lizard.predict[,3]

summary.df <- unique(df.2[,-3])
names(summary.df)[1:3] <- c('Sex','Location','Means')


qplot(data=summary.df,
        y=Means,
        x=Location,
        group=Sex,
        ymin=lower,
        ymax=upper,
        geom=c("point", "errorbar", "line"),
        color=Sex,
        shape=Sex,
        width=0.25)   theme_bw()
  

PS: Если я кажусь резким здесь и там, это не предназначено. Английский не мой родной язык, и я до сих пор не знаком с тонкостями языка.

[Потенциальное бесстыдное продвижение] Вам следует рассмотреть функции compareCats и rxnNorm в пакете HandyStuff, доступном по адресу www.github.com/bryanhanson/HandyStuff Предупреждение: я не уверен, что это работает без проблем с R 2.14. В частности, rxnNorm выглядит как график, который вы пытаетесь создать, плюс он дает вам множество вариантов с точки зрения обобщающей статистики и оформления графика. Но для этого необходимо, чтобы ваши ученики установили отдельный пакет, поэтому, возможно, вы его исключите (но это позволяет учащимся сосредоточиться на представлении и анализе данных). График из приведенного здесь примера ?rxnNorm.

С rxnNorm у вас есть выбор из нескольких способов вычисления CI, управляемых аргументом «метод». Вот фактические функции (из пакета ChemoSpec).

 > seX <- function (x)  sd(x, na.rm = TRUE)/sqrt(length(na.omit(x)))
> <environment: namespace:ChemoSpec>
> 
> seXy <- function (x)  {
>     m <- mean(na.omit(x))
>     se <- seX(x)
>     u <- m   se
>     l <- m - se
>     c(y = m, ymin = l, ymax = u) } <environment: namespace:ChemoSpec>
> 
> 
> seXy95 <- function (x)  {
>     m <- mean(na.omit(x))
>     se <- seX(x)
>     u <- m   1.96 * se
>     l <- m - 1.96 * se
>     c(y = m, ymin = l, ymax = u) } <environment: namespace:ChemoSpec>
> 
> 
> seXyIqr <- function (x)  {
>     i <- fivenum(x)
>     c(y = i[3], ymin = i[2], ymax = i[4]) } <environment: namespace:ChemoSpec>
> 
> seXyMad <- function (x)  {
>     m <- median(na.omit(x))
>     d <- mad(na.omit(x))
>     u <- m   d
>     l <- m - d
>     c(y = m, ymin = l, ymax = u) } <environment: namespace:ChemoSpec>