Новичок пробует модель линейных смешанных эффектов в студии R - ПОЛНАЯ НЕУДАЧА

После более чем часа поиска (этот форум, Youtube, заметки о занятиях, Google) я не нашел помощи по моему вопросу. Я полный новичок, который ничего не знает ни о R, ни о статистике.

Я пытаюсь создать линейную модель смешанных эффектов в R. Я измеряю ширину листа в трех разных местах (Джексонвилл, Флорида, Огаста, Джорджия, и Атланта, Джорджия), и в этих трех местах есть растения с высоким и низким содержанием азота. азотный участок. У меня есть 150 измерений листьев с 50 деревьев.

Мое ограниченное понимание говорит мне, что ширина листа — это непрерывная переменная отклика, а город и участок — дискретные независимые переменные. Случайным эффектом будут отдельные деревья, поскольку ширина листа в пределах одного дерева не является независимой.

Я использовал "nlme" для создания модели:

leaf.width.model <- lme(width ~ city*plot, (1|tree.id), data=leaf)

Затем я провел тест ANOVA, и он предположил, что что-то происходит с городом и взаимодействием между городом и участком. Вот где я застрял. Я хочу сделать сюжет, в котором есть линии для всех трех городов, но я понятия не имею, как это сделать. Когда я пытаюсь использовать функцию сюжета, я просто получаю коробочную диаграмму.

Я буквально часами пытался, и я еще больше потерян и сбит с толку, чем раньше.

1) Как я могу сделать этот график?

2) Какие еще тесты я должен сделать, чтобы проанализировать и/или визуализировать эти данные?

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

Спасибо,

Богатый

P.S Вот вывод функции dput:

> dput(tree) structure(list(tree.id = structure(c(24L, 24L, 32L, 25L, 25L, 24L, 24L, 32L, 25L, 25L, 43L, 45L, 45L, 43L, 23L, 23L, 45L, 45L, 23L, 23L, 41L, 41L, 38L, 11L, 11L, 38L, 41L, 41L, 11L, 11L, 14L, 14L, 29L, 13L, 13L, 14L, 14L, 29L, 13L, 13L, 4L, 4L, 1L, 1L, 20L, 1L, 1L, 20L, 6L, 8L, 8L, 5L, 5L, 6L, 4L, 4L, 8L, 8L, 5L, 5L, 9L, 9L, 10L, 10L, 12L, 12L, 13L, 13L, 22L, 22L, 23L, 23L, 24L, 24L, 25L, 25L, 25L, 25L, 40L, 40L, 41L, 41L, 38L, 38L, 39L, 39L, 14L, 14L, 14L, 15L, 15L, 28L, 28L, 29L, 29L, 35L, 35L, 36L, 36L, 37L, 37L, 42L, 42L, 43L, 43L, 44L, 44L, 45L, 45L, 46L, 46L, 47L, 47L, 2L, 1L, 3L, 3L, 4L, 4L, 7L, 11L, 11L, 16L, 16L, 20L, 20L, 21L, 21L, 17L, 17L, 18L, 18L, 19L, 19L, 26L, 26L, 27L, 27L, 30L, 30L, 31L, 31L, 32L, 32L, 33L, 33L, 34L, 34L, 48L), .Label = c("Tree_112", "Tree_112 ", "Tree_115", "Tree_130", "Tree_137", "Tree_139", "Tree_140", "Tree_141", "Tree_153", "Tree_154", "Tree_156", "Tree_159", "Tree_166", "Tree_169", "Tree_171", "Tree_180", "Tree_182", "Tree_184", "Tree_185", "Tree_202", "Tree_213", "Tree_218", "Tree_222", "Tree_227", "Tree_239", "Tree_242", "Tree_246", "Tree_247", "Tree_252", "Tree_260", "Tree_267", "Tree_269", "Tree_271", "Tree_272", "Tree_291", "Tree_293", "Tree_298", "Tree_327", "Tree_329", "Tree_336", "Tree_350", "Tree_401", "Tree_403", "Tree_405", "Tree_407", "Tree_409", "Tree_420", "Tree_851"), class = "factor"), city = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 2L, 2L, 3L, 3L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("Atlanta", "Augusta", "Jacksonville"), class = "factor"),     plot = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("High-N", "Low-N"), class = "factor"), width = c(0.66, 0.716, 0.682, 0.645, 0.645, 0.696, 0.733,
0.707, 0.668, 0.686, 0.617, 0.733, 0.73, 0.615, 0.669, 0.746, 0.687, 0.682, 0.76, 0.713, 0.651, 0.664, 0.679, 0.729, 0.756, 
0.669, 0.647, 0.713, 0.767, 0.685, 0.69, 0.731, 0.781, 0.729, 
0.725, 0.739, 0.769, 0.791, 0.676, 0.688, 0.719, 0.753, 0.748, 
0.791, 0.785, 0.78, 0.723, 0.756, 0.664, 0.645, 0.653, 0.615, 
0.591, 0.642, 0.693, 0.716, 0.694, 0.676, 0.662, 0.629, 0.665, 
0.748, 0.726, 0.693, 0.715, 0.714, 0.764, 0.732, 0.61, 0.721, 
0.703, 0.713, 0.746, 0.752, 0.662, 0.733, 0.707, 0.674, 0.734, 
0.79, 0.732, 0.794, 0.703, 0.712, 0.737, 0.731, 0.747, 0.746, 
0.787, 0.709, 0.716, 0.764, 0.77, 0.764, 0.802, 0.663, 0.777, 
0.642, 0.779, 0.81, 0.724, 0.645, 0.68, 0.637, 0.695, 0.768, 
0.761, 0.7, 0.759, 0.726, 0.696, 0.794, 0.774, 0.799, 0.747, 
0.606, 0.691, 0.733, 0.707, 0.698, 0.706, 0.72, 0.694, 0.697, 
0.737, 0.716, 0.73, 0.706, 0.667, 0.734, 0.528, 0.695, 0.684, 
0.763, 0.733, 0.809, 0.6, 0.676, 0.718, 0.759, 0.609, 0.665, 
0.667, 0.647, 0.701, 0.663, 0.688, 0.693, 0.899)), .Names = c("tree.id", "city", "plot", "width"), class = "data.frame", row.names = c(NA, -149L))

Большое спасибо всем за ваши комментарии, я искренне ценю помощь каждого!


r mixed-models random-effects
person Richard Gourderton    schedule 22.03.2018    source источник
comment
Привет, Ричард, добро пожаловать в SO! Можете ли вы предоставить нам снимок ваших данных, используя, например, dput? Это облегчит помощь (подробнее см. здесь: stackoverflow.com/questions/5963269/). Кроме того, вам не нужно быть самоуничижительным, чтобы получить здесь помощь. На самом деле, все, что нас волнует, — это содержание вашего вопроса, а не то, как он возник (если только это не актуально), и не то, насколько вы хороши в R или статистике.   -  person tifu    schedule 22.03.2018
comment
Я не уверен, что линейный сюжет имеет смысл. Размещение непрерывного отклика на оси Y понятно, но у вас нет ничего непрерывного, что можно было бы разместить на оси X. Мне также неясна ваша цель для графика: вы хотите построить необработанные данные или прогноз модели (или оба)? Если это необработанные данные, может иметь смысл блочная диаграмма, возможно, с разными цветами боксов для каждого графика.   -  person Gregor Thomas    schedule 22.03.2018

Ответы (2)


arrow_upward
1
arrow_downward

Как указано в комментариях, линейный график может не иметь смысла для ваших данных, поскольку вы изучаете, как ширина изменяется в отдельных категориях (в отдельных городах и на отдельных участках). Боксплоты имеют смысл, поскольку вы можете создавать их для каждого из взаимодействий города и сюжета. Чтобы дать вам представление о том, что вы можете сделать, я сгенерировал некоторые фальшивые данные и сделал пример графика, который может быть вам полезен:

# fake data
leaf <- data.frame(tree.id = rep(1:50, each = 3), 
                   city = rep(c("Jackson", "Augusta", "Atlanta"), each = 50),
                   plot = rep(1:6, each = 25))
# I'll make the average of width different for each plot
leaf$width <- rnorm(nrow(leaf), leaf$plot, 1)

# plotting the data
library(ggplot2) # this is a great library for plotting in R

ggplot(leaf, aes(x = factor(plot), y = width, color = factor(plot))) + 
  facet_grid(~city, scales = 'free_x') + # This creates a subplot for each city
  geom_boxplot() + 
  geom_point(position = "jitter") + 
  theme_bw()

На этом графике я добавил точки (ширину листа для каждого отдельного дерева), но я «поколебал» их, то есть слегка нарушил их положение, чтобы они не накладывались друг на друга и были видны все. Вы можете удалить это, если хотите.

Исследовательский анализ данных должен приносить удовольствие! И я думаю, что визуализация — это хорошее начало, когда вы начинаете заниматься статистикой. Надеюсь, это окажется полезным для вас.

person gfgm    schedule 22.03.2018
comment
gfgm, большое спасибо за помощь! Есть ли способ представить среднее значение каждой группы чем-то вроде треугольника или звезды? Я искал, как добавить среднее значение в ggplot, но не смог найти ничего, что бы сработало. Еще раз спасибо. - person Richard Gourderton; 30.03.2018
comment
вы можете поместить большой черный крест в среднее значение каждой группы, добавив + stat_summary(fun.y = "mean", geom = "point", color = 'black', size=10, shape = "+") в конец вызова функции построения графика. Если ответ полезен, проголосуйте за него! - person gfgm; 30.03.2018
comment
Последняя вещь! Как мне изменить порядок диаграмм, чтобы они переместились из Джексонвилля в Огасту и в Атланту? Я попытался использовать функцию конкатенации c, но она либо просто удаляет все города, либо создает отдельный объект переменной с именем city, в котором есть только 3 элемента (только города). Большое вам спасибо! Я проголосовал и отметил как лучший ответ - person Richard Gourderton; 05.04.2018

arrow_upward
0
arrow_downward 

leaf.width.model <- lme(width ~ city*plot, (1|tree.id), data=leaf)

В этой модели, если вы хотите что-то нарисовать, вы, вероятно, пытаетесь ответить: сколько составляет средняя ширина листа для всех деревьев в каждом городе для каждого типа участка.

Чтобы показать эту информацию на рисунке, вам нужно отложить width по оси y, построить plot(высокий и низкий уровень азота) по оси x и сгруппировать данные по city. Тогда вы получите 3 строки, о которых вы говорите. Однако вам нужно получить среднюю ширину в каждой группе, так как вы хотите показать только городские различия.

Чтобы получить этот график из необработанных данных: (используя поддельные данные, предоставленные gfgm)

set.seed(100)
leaf <- data.frame(tree.id = rep(1:50, each = 3), 
                   city = rep(c("Jackson", "Augusta", "Atlanta"), each = 50),
                   plot = rep(c(1, 0), each = 25))

# I'll make the average of width different for each plot
leaf$width <- rnorm(nrow(leaf), leaf$plot, 1)

library(plotly)
library(tidyverse)
leaf %>% 
  group_by(city,plot) %>% 
  summarise(avwidth = mean(width, na.rm=T),
            avsd = 1.96*sd(width, na.rm=T)/sqrt(25)) %>%
  plot_ly(x = ~plot, y = ~avwidth, color= ~city, 
          type="scatter", mode="markers+lines", 
          error_y = ~list(array=avsd)
          )

person Rajesh Talluri    schedule 22.03.2018
comment
Я только что попробовал это, и это работает потрясающе! Итак, насколько я понимаю, точки представляют собой среднюю ширину для этой группы (например, Атланта, High-N). Есть ли способ добавить полосу ошибок к точкам? Большое спасибо Раджеш Таллури! - person Richard Gourderton; 30.03.2018
comment
Я отредактировал ответ, чтобы добавить стандартные ошибки. Здесь я просто использую обычные симметричные стандартные ошибки среднего значения +-1,96 * sd / sqrt (n) - person Rajesh Talluri; 03.04.2018
comment
Большое спасибо! Я проголосовал за ваш ответ (хотя он говорит, что он не будет отображаться, потому что я новичок). - person Richard Gourderton; 05.04.2018