Новичок, пытающийся моделировать линейные смешанные эффекты в студии R - TOTAL FAIL

Question

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

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

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

Я использовал "NLME", чтобы сделать модель:

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

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

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

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

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

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

Спасибо,

Богатые

PS Вот выход из 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))

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

Answer 1

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

# 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()

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

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

Answer 2

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

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

Чтобы показать эту информацию на рисунке, вам нужно построить width на графике оси Y plot(высокий и низкий азот) на оси х и сгруппировать данные по 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)
          )