Интерпретация результатов анализа данных

Question

Мне нужна помощь в интерпретации результатов моих данных для моей дипломной работы. Я изучаю влияние лесопользования на размножение певчих птиц. Я взял данные, которые были собраны за последние шесть лет, и провел анализ AI Cc с использованием моделей множественной линейной регрессии со случайными и фиксированными эффектами. Я получил значения AI Cc для каждой модели и определил, какие модели лучше всего описывают различия в данных, но у меня возникают проблемы с переводом этих результатов в представительный формат с использованием графиков, рисунков и, ну, в общем, слов. Как мне показать, каковы результаты моего анализа в сжатой форме, и есть ли способы составить график такого анализа, чтобы он хорошо выглядел в статье? Я открыт для использования R для отображения своих результатов, но я все еще очень плохо знаком с технологией, и поэтому я не знаю, какие варианты есть в R, или какие шаги мне нужно предпринять, чтобы сделать графики презентабельными.

Answer 1

Использование библиотеки ggplot2 может создавать некоторые очень публикуемые графики непосредственно из выходных данных LMM и GLMM из объектов lme4 :: lmer и glmer, но имеет кривую обучения. Если у вас есть как непрерывные, так и категориальные предикторы, вот что я использовал для текущего проекта. Там много настроек, которые выглядят пугающе, но это вопрос поиска снова и снова решения и поиска правильного куска кода, который решает вашу проблему. Если у вас есть вопросы, дайте мне знать. Я немного упомянул, чтобы помочь.

Я хотел бы знать, изменилось ли значение алкоголя hedoni c для крыс, выращенных в разных средах (Среда) и при разных уровнях воздействия алкоголя в подростковом возрасте (Условие ).

У меня было 4 переменные "c .con c" (средняя концентрация, центрированная во избежание проблем с инфляцией дисперсии из-за мультиколинейности), "Пол", "Состояние", "Окружающая среда". Моя переменная концентрации была в пределах субъекта и, следовательно, моя повторная мера.

#Load in the required libraries

library("MASS")
library("lattice")
library("boot")
library("car")
library("emmeans")
library("lme4")
library("zoo")
library("tidyr")
library("multcomp")
library("foreign")
library("msm")
library("ggplot2")
library("effects")
library("lmerTest")

#Run the model and put the results into an object i called "Ehed" for Ethanol Hedonics

Ehed <-glmer(Total.Hedonic ~ c.conc*Sex*Condition*Environment
                 + (c.conc|RatID), data=mydetoh, family=poisson)
    summary(Ehed)

#Always check the normality of your residuals so you don't violate assumptions of residual distributions.
#In my case, they were very normal and the graph below will be saved as a .png to my current working directory to visually compare my residual distribution to a normal curve.
#to view your current working directory enter "getwd()" 

    #Residual Graph
      #make PNG file
      png("COBRE-2 Ehed Res Plot.png", width = 300, height = 300)
      #plot residual density function
      plot(density(residuals(Ehed)), 
           main="", xlab="", frame= FALSE)
      #Add normal distribution to the residual plot for comparison
      Ehed.res = residuals(Ehed)
      Ehed.m = mean(Ehed.res)
      Ehed.std = sqrt(var(Ehed.res))
      curve(dnorm(x, mean=Ehed.m, sd=Ehed.std), col="darkblue", lwd=2, add=TRUE, yaxt="n")
      #close the file
      dev.off()

После того, как вы запустите свои модели, вам нужно будет сделать несколько вещей, чтобы убедиться, что вы не столкнетесь с проблемами. Настройте свои разрывы оси на основе ваших максимальных значений вашей прогнозируемой переменной. Измените масштаб и отцентрируйте все непрерывные переменные, которые должны были быть отцентрированы перед анализом, чтобы избежать проблем с мультиколинейностью (это делается на том же шаге, что и следующая вещь в этом списке). Вытащите свои эффекты из объекта модели, чтобы эффективно рассчитать ленты ошибок для стандартной ошибки среднего (SEM).

#Predicted Graphs####

##Graph Setup####

##ETHANOL####
  #Axis and Break/Label Setup
    #Y Axis Breaks/Labels
      # generate hedonic (my predicted variable) Y axis break positions
      Ehed.ybreaks = c(0,50,100,150,200,250,300,350,400)
      # and Y labels; the ,"",100... omits the label for 50 but leaves the tick mark.
      # The length of the vectors must be the same so the ""s are necessary for this trick.
      Ehed.ylabels = as.character(c(0,"",100,"",200,"",300,"",400))

    #X Axis Breaks/Labels
      #assign X break positions for Concentration to object. My Concentration variable was 5%, 10%, 20%, etc... and appears below in the list (c())
      E.xbreaks = c(5,10,20,30,40)
      #assigns the values of the breaks to the X labels
      E.xlabels = as.character(E.xbreaks)

    ###Ethanol Hedonic Graph______________________________________________________
      #pull the effects from the GLMER model object & calculate confidence intervals for graphing
      Ehed.eff <- Effect(c("c.conc","Sex","Condition","Environment"),Ehed,
                         #se is std err and the level is the confidence level. .68 = actual std err for conf int. lower and upper.
                         se=list(level=.68),
                         #the xlevels command is used to increase the number of points calculated to smooth the error ribbons to look more curved (default = 5). I also center my concentration variable here to line up with the numbers that were analyzed as centered variables and remove this in the next step so everything is at their original values.
                         xlevels=list(c.conc=c(.05-mean.etoh.conc,
                                               .075-mean.etoh.conc,
                                               .10-mean.etoh.conc,
                                               .125-mean.etoh.conc,
                                               .15-mean.etoh.conc,
                                               .175-mean.etoh.conc,
                                               .20-mean.etoh.conc,
                                               .225-mean.etoh.conc,
                                               .25-mean.etoh.conc,
                                               .275-mean.etoh.conc,
                                               .30-mean.etoh.conc,
                                               .325-mean.etoh.conc,
                                               .35-mean.etoh.conc,
                                               .375-mean.etoh.conc,
                                               .40-mean.etoh.conc)))
#Writing Ehed.eff to a data frame to more easily use it with other functions later.
      Ehed.eff.df <-as.data.frame(Ehed.eff)
      #Converting back from the centering and rescaling.
      Ehed.eff.df$Concentration <- (Ehed.eff.df$c.conc+mean.etoh.conc)*100
      #Instead of trying to relabel these values in the ggplot2 object just recode them here and save some trouble
      Ehed.eff.df$Sex <-car::recode(Ehed.eff.df$Sex, "'F' = 'Female'; 'M' = 'Male'")
#The mydetoh is the original data set. If i want to show points for each individual I can use this data set layered on top of my other graph.
      mydetoh$Sex <-car::recode(mydetoh$Sex, "'F' = 'Female'; 'M' = 'Male'")

      #Check that everything looks right.
      View(Ehed.eff.df)
      View(mydetoh)

      #make a new file
      png("Fig Sample Ethanol Hedonic lines.png", width = 800, height = 600)
      #Start your plot and write it to an object for later reference. I chose Ehed.ggp because it is the Ehed model's ggplot. The fit variable below is generated when you pull the effects from the model. It is the predicted values.
#Because I am spliting the plot into two panes by sex, only Condition and Environment need to appear in my group, col (color), fill, and linetype arguments. The names must match the names of your variables from your model EXACTLY. The scale_color_manual etc... all align with these values. I used Hex colors, you can also just type "red" with the quotes instead.

      Ehed.ggp <-ggplot(Ehed.eff.df,
                        aes(Concentration,fit,
                            group=interaction(Condition,Environment),
                            col=interaction(Condition,Environment),
                            fill=interaction(Condition,Environment),
                            linetype=interaction(Condition,Environment),
                            shape=interaction(Condition,Environment)))+
        #adds each individual's points to the data. Leave this commented out if you dont need to do that.
        #geom_point(data=mydetoh,aes(x=Concentration, y=Total.Hedonic),stroke=1.5,size=4,alpha=0.60)+
        geom_smooth(data=Ehed.eff.df, se=FALSE, method="glm", method.args = list(family = "poisson"),size=1.5)+
        ## colour=NA suppresses edges of the ribbon
        geom_ribbon(data=Ehed.eff.df,colour=NA,alpha=0.25,
                    aes(ymin=lower,ymax=upper))+
        #labs(tag="A.")+  #If you do not need a panel tag (e.g. A., B., C. etc...) for a graph that will become part of a larger plot, comment this command out
        scale_color_manual("",values=c("#0000ff", "#7d7dff","#ff0000","#ff7d7d","#000000","#808080"), labels=c('EC+ETOH','EC+SAL','IC+ETOH','IC+SAL','SC+ETOH','SC+SAL'))+
        scale_fill_manual("",values=c("#0000ff", "#7d7dff","#ff0000","#ff7d7d","#000000","#808080"), labels=c('EC+ETOH','EC+SAL','IC+ETOH','IC+SAL','SC+ETOH','SC+SAL'))+
        scale_linetype_manual("",values=c("solid","twodash","solid","twodash","solid","twodash"), labels=c('EC+ETOH','EC+SAL','IC+ETOH','IC+SAL','SC+ETOH','SC+SAL'))+
        scale_shape_manual("",values=c(15,0,16,1,17,2), labels=c('EC+ETOH','EC+SAL','IC+ETOH','IC+SAL','SC+ETOH','SC+SAL'))+
        scale_x_continuous(expand=c(0,0), limits = c(0,42), breaks=E.xbreaks, labels=E.xlabels)+
        scale_y_continuous(expand=c(0,0), limits = c(0,410), breaks=Ehed.ybreaks, labels=Ehed.ylabels)+
        theme_classic()+
        theme(strip.background = element_rect(colour="white"),
              strip.text.x = element_text(size=18,face="bold"),
              panel.spacing = unit(1,"cm"),
              axis.title = element_text(size=22),
              axis.text = element_text(size=21,color="black",face="bold"),
              axis.line = element_line(size=1.3),
              axis.ticks = element_line(size=1.3, color="black"),
              axis.ticks.length = unit(0.2, "cm"),
              axis.title.y = element_text(margin = margin(t = 0, r = 18, b = 0, l = 0)),
              axis.title.x = element_text(margin = margin(t = 13, r = 0, b = 0, l = 0)),
              legend.title = element_blank(),
              legend.text = element_text(size=18, face="bold"),
              legend.justification = "top",
              legend.key.size = unit(1, "cm"),
              #legend.position = c(0.75, .85),
              #plot.tag = element_text(size=36, face="bold"),
              plot.tag.position = c(0.05, 0.95))+
        facet_grid(. ~ Sex)+
        xlab("Ethanol % (v/v)")+
        ylab("Hedonic Responses (+/-SEM)")

      Ehed.ggp

      #close the file
      dev.off()

К сожалению, у меня нет времени, чтобы объяснить больше прямо сейчас. Надеюсь, это поможет.