解决办法：
首先导入&&子集化数据
#I called mine Ancovas. --> Note, export your df as .csv to work with it in R.
Ancovas <- read.csv("~/Dropbox/YOUR DATAFILE NAME.csv")
#Next, subset your data by the two conditions (e.g. "l"=light, "d"=dark), and both treatments (e.g. "MQ"=water, "DOM"=media)
AncovasL <- Ancovas[(Ancovas$UV == "Light"), ]
AncovasL.MQ <- AncovasL[(AncovasL$DOM == "MQ"), ]
AncovasL.DOM <- AncovasL[(AncovasL$DOM == "DOM"), ]
AncovasD <- Ancovas[(Ancovas$UV == "Dark"), ]
AncovasD.MQ <- AncovasD[(AncovasD$DOM == "MQ"), ] #This code only keeps what is inside the brackets
AncovasD.DOM <- AncovasD[(AncovasD$DOM == "DOM"), ] #Note, adding and "!" after square bracket removes what is in " ".
创建回归函数
#--> this code was gathered from several sites.
#note: I don't understand the logic of how the numbers in brackets are organized. But this essentially pulls some information from the fit model. i.e. [9] means find the 9th value in the list (I think)
regression = function(Ancovas){
fit <- lm(AvgBio ~ Exposure, data=Ancovas)
slope <- round(coef(fit)[2],1)
intercept <- round(coef(fit)[1],0)
R2 <- round(as.numeric(summary(fit)[8]),3)
R2.Adj <- round(as.numeric(summary(fit)[9]),3)
p.val <- signif(summary(fit)$coef[2,4], 3)
c(slope,intercept,R2,R2.Adj, p.val) }
现在通过TREATMENT拆分回归数据并应用回归函数
#Call your column "Treatments"
regressions_dataL.MQ <- ddply(AncovasL.MQ, "Treatment", regression) #For light samples using water
regressions_dataL.DOM <- ddply(AncovasL.DOM, "Treatment", regression) #For light samples using media
regressions_dataD.MQ <- ddply(AncovasD.MQ, "Treatment", regression) #For dark samples using water
regressions_dataD.DOM <- ddply(AncovasD.DOM, "Treatment", regression) #For dark samples using media
#Rename columns
colnames(regressions_dataL.MQ) <-c ("Treatment","slope","intercept","R2","R2.Adj","p.val")
colnames(regressions_dataL.DOM) <-c ("Treatment","slope","intercept","R2","R2.Adj","p.val")
colnames(regressions_dataD.MQ) <-c ("Treatment","slope","intercept","R2","R2.Adj","p.val")
colnames(regressions_dataD.DOM) <-c ("Treatment","slope","intercept","R2","R2.Adj","p.val")
为数字创建主题
#Yes I like to hyper control every aspect of my theme
theme_new <- theme(panel.background = element_rect(fill = "white", linetype = "solid", colour = "black"),
legend.key = element_rect(fill = "white"), panel.grid.minor = element_blank(), panel.grid.major = element_blank(),
axis.text.x=element_text(size = 11, angle = 0, hjust=0.5), #axis numbers (set it to 1 to place it on left side, 0.5 for middle and 0 for right side)
axis.text.y=element_text(size = 13, angle = 0),
plot.title=element_text(size=15, vjust=0, hjust=0), #hjust 0.5 to center title
axis.title.x=element_text(size=14), #X-axis title
axis.title.y=element_text(size=14, vjust=1.5), #Y-axis title
legend.position = "top",
legend.title = element_text(size = 11, colour = "black"), #Legend title
legend.text = element_text(size = 8, colour = "black", angle = 0), #Legend text
strip.text.x = element_text(size = 9, colour = "black", angle = 0), #Facet x text size
strip.text.y = element_text(size = 9, colour = "black", angle = 270)) #Facet y text size
guides_new <- guides(color = guide_legend(reverse=F), fill = guide_legend(reverse=F)) #Controls the order of your legend
Colours <-
rainbow_hcl(length(levels(factor(StackedTable$DOM))), start = 30, end = 300) #Yes I am Canadian so Colours has a "u"
Colours[5] <- "#47984c" #Green
Colours[4] <- "#7b64b4" #Purple-grey
Colours[3] <- "#ff7f50" #Orange
Colours[2] <- "#cc3636" #Red
Colours[1] <- "#4783ba" #Blue
创建两个稍后合并的数字
#Making plot for panel A ("Dark condition")
PlotA <-
ggplot(AncovasD, aes(x=as.numeric(Time.h), y=as.numeric(Measurement), fill=as.factor(Treatment))) +
geom_smooth(data=subset(AncovasD,Treatment =="MQ"), aes(Time.h,Measurement,color=factor(Treatment)),method="lm", formula = y~x, se=T, show.legend = F) +
geom_smooth(data=subset(AncovasD,Treatment =="DOM"), aes(Time.h,Measurement,color=factor(Treatment)),method="lm", formula = y~x, se=T, show.legend = F) + #You need this line twice, once for each condition
geom_errorbar(data=AncovasD, aes(ymin=Measurement-SD, ymax=Measurement+SD), width=0.2, colour="#73777a", size = 0.5) + #Change width based on the size of your X-axis
geom_point(shape = 21, size = 3, colour = "black", stroke = 1) + #colour is the outline of the circle, stroke is the thickness of that outline
facet_grid(Treatment ~ UV) + #This places all your treatments into a grid. Change the order if you want them horizontal. Use "." if you do not want a label.
geom_label(data=regressions_dataD.MQ, inherit.aes=FALSE, size=0.7, colour=Colours[1], #Add label for DOM regressions, specify same colour as your legend, change size depending on how large you want the text
aes(x=-0.1, y=41, label=paste(" ", "m == ", slope, "\n " , #replace this line with the values you want: e.g. R-squared=("R2 == ", R2.Adj) ; intercept=("b == ", intercept). The "\n " makes a second line
" ", "p == ", p.val ))) + #This completes the first label. Repeat same process for second label.
geom_label(data=regressions_dataD.DOM, inherit.aes=FALSE, size=0.7, colour=Colours[2],
aes(x=-0.1, y=4, label=paste(" ", "m == ", slope, "\n " ,
" ", "p == ", p.val )))
#Now for the irradiated samples "light" plot (Panel B)
PlotB <-
ggplot(AncovasL, aes(x=as.numeric(Time.h), y=as.numeric(Measurement), fill=as.factor(Treatment))) + #Same as above but use your second dataframe.
geom_smooth(data=subset(AncovasL,Treatment =="MQ"), aes(Time.h,Measurement,color=factor(Treatment)),method="lm", formula = y~x, se=T, show.legend = F) +
geom_smooth(data=subset(AncovasL,Treatment =="DOM"), aes(Time.h,Measurement,color=factor(Treatment)),method="lm", formula = y~x, se=T, show.legend = F) +
geom_errorbar(data=AncovasL, aes(ymin=Measurement-SD, ymax=Measurement+SD), width=0.2, colour="#73777a", size = 0.5) +
geom_point(shape = 21, size = 3, colour = "black", stroke = 1) +
facet_grid(Treatment ~ UV) +
geom_label(data=regressions_dataL.MQ, inherit.aes=FALSE, size=0.7, colour=Colours[1],
aes(x=-0.1, y=41, label=paste(" ", "m == ", slope, "\n " ,
" ", "p == ", p.val ))) +
geom_label(data=regressions_dataL.DOM, inherit.aes=FALSE, size=0.7, colour=Colours[2],
aes(x=-0.1, y=4, label=paste(" ", "m == ", slope, "\n " ,
" ", "p == ", p.val )))
