解读数据分析的结果

Question

我正在寻找一些帮助解释我的数据的结果，为我的研究生论文。我正在研究森林管理对歌鸟繁殖的影响。我采用了过去六年收集的数据，并使用具有随机和固定影响的多个线性回归模型进行了AICc分析。我已经获得了每个模型的AICc值，并确定了哪些模型最能描述数据中的变化，但我无法使用图形、图形和文字将这些结果转换为可呈现的格式。我如何以简明的方式显示我的分析结果，有没有方法用图表表示这种类型的分析，以便在论文中看起来很好？我对使用R绘制结果持开放态度，但我对这项技术仍然非常陌生，因此我不知道R中有哪些选项，也不知道我需要采取哪些步骤才能使绘图更具表现力。

Answer 1

使用ggplot2库可以直接从lme4::lmer和glmer对象的LMM和GLMM的输出中生成一些非常可发布的图形，但确实有一个学习曲线。如果你同时有连续预测和分类预测，这里是我在当前项目中使用的。有许多设置看起来令人生畏，但这是一个反复搜索解决方案并找到解决问题的正确代码的问题。如果你有什么问题，请告诉我。我已经做了相当多的注释来提供帮助。

我想知道在不同环境(环境)和青春期不同程度酒精暴露的大鼠中，酒精的享乐性价值是否发生了变化。

我有4个变量"c.conc“(集中均值，以避免由于多重共线性而导致的方差膨胀问题)，”性别“，”条件“，”环境“。我的注意力变量在受试者内部，因此，我的重复测量。

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

不幸的是，我现在没有时间解释更多。希望这能有所帮助。