我正在寻求一些帮助来解释我的研究生论文的数据结果。我正在研究森林管理对鸣禽繁殖的影响。我采用了过去六年收集的数据，并使用具有随机和固定效应的多重线性回归模型进行了 AICc 分析。我已经获得了每个模型的 AICc 值，并确定了哪些模型最能描述数据的变化，但是我无法使用图形、数字和文字将这些结果转化为可呈现的格式。我如何以简洁的方式显示我的分析结果，有没有办法将这种类型的分析绘制成图表，以便在论文中看起来不错？我愿意使用 R 来绘制我的结果，但我对这项技术仍然很陌生，所以我不知道 R 中有哪些选项，或者我需要采取哪些步骤才能使图表看起来像样。
【问题讨论】：
标签： plot
使用 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()
很遗憾，我现在没有时间解释更多。希望这会有所帮助。
【讨论】：