All Morphology: Wing LOESS
Anastasia Bernat
4/8/2021
#knitr::purl("wing_loess.Rmd", output = "wing_loess.R") # convert Rmd to R scriptCleaning Data and Sourcing Scripts
source_path = "~/Desktop/git_repositories/SBB-dispersal/avbernat/Rsrc/"
script_names = c("compare_models.R",
"regression_output.R",
"clean_morph_data.R",
"AICprobabilities.R",
"remove_torn_wings.R",
"get_Akaike_weights.R")
for (script in script_names) {
path = paste0(source_path, script)
source(path)
}##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric##
## Attaching package: 'lubridate'## The following object is masked from 'package:cowplot':
##
## stamp## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, union##
## Attaching package: 'dplyr'## The following object is masked from 'package:gridExtra':
##
## combine## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionRead the data
data_list <- read_morph_data("data/allmorphology05.10.21-coors.csv")## morph types: L S NA LS L/S SL
## recoding missing morph types...
## S if wing2thorax <=2.2, L if wing2thorax >=2.5
##
## ambiguous wing morph bug count: 32raw_data = data_list[[1]]
data_long = data_list[[2]]
all_bugs = nrow(data_long)
data_long = remove_torn_wings(data_long)##
## number of bugs with torn wings: 193Wing2Body
SE = function(x){sd(x)/sqrt(length(x))}
w2b_summary<-aggregate(wing2body~sex*pophost*dates, data=data_long, FUN=mean)
w2b_summary$se<-aggregate(wing2body~sex*pophost*dates, data=data_long,
FUN=SE)$wing2body
#jitter slightly
jitter = runif(n=nrow(w2b_summary), min=-0.5, max=0.5)
w2b_summary$dates <- w2b_summary$dates + jitter# check for residuals for loess
l1 = lowess(w2b_summary$dates, w2b_summary$wing2body, f=0.4)
plot(w2b_summary$dates, w2b_summary$wing2body)
lines(l1, type = "l")
plot_lowess_residuals = function(lfit, x, y) {
lfun <- approxfun(lfit)
fitted <- lfun(x)
resid <- y-fitted
plot(fitted,resid)
abline(h=0,col=8)
}
plot_lowess_residuals(l1, w2b_summary$dates, w2b_summary$wing2body)
d = w2b_summary
# general plotting features
xlab_years = na.omit(sort(unique(data_long$dates))[-2])
host_colors_shade = c("turquoise3", "green")
host_colors_pts = c("turquoise3", "springgreen4")
# date plotting features
threshold = unique(d$dates)[13]
threshold## [1] "2016-12-01"d$date_b = "2013-2014"
d$date_b[d$dates >= threshold] = "2015-2020"
d$date_b = as.factor(d$date_b)
w2b_table<-aggregate(wing2body~dates, data=data_long, FUN=mean)
w2b_table$date_b = "2013-2014"
w2b_table$date_b[w2b_table$dates >= threshold] = "2015-2020"
w2b_table$date_b = as.factor(w2b_table$date_b)
last_n_years = w2b_table[w2b_table$date_b=="2015-2020",]
last_n_years## dates wing2body date_b
## 5 2017-08-01 0.7331342 2015-2020
## 6 2018-09-01 0.7301813 2015-2020
## 7 2019-05-01 0.7341616 2015-2020
## 8 2019-10-01 0.7299658 2015-2020
## 9 2020-02-01 0.7234360 2015-2020dfF = d[d$sex=="F",]
dfM = d[d$sex=="M",]
dfF$pophost[dfF$pophost=="C.corindum"]<-"C. corindum"
dfF$pophost[dfF$pophost=="K.elegans"]<-"K. elegans"
dfM$pophost[dfM$pophost=="C.corindum"]<-"C. corindum"
dfM$pophost[dfM$pophost=="K.elegans"]<-"K. elegans"
dfF$pophost = factor(dfF$pophost, levels = c("K. elegans", "C. corindum") )
dfF$`Host Plant` = dfF$pophost
dfM$pophost = factor(dfM$pophost, levels = c("K. elegans", "C. corindum") )
dfM$`Host Plant` = dfM$pophostp0 = ggplot() +
ggtitle("C") + xlab("Year") + ylab("Wing-to-Body Ratio") +
geom_vline(xintercept = xlab_years, color="gainsboro") +
geom_smooth(data=w2b_table, method="lm", se=FALSE, linetype = "dashed",
mapping = aes(x = dates, y = wing2body), colour="black", lwd=0.5) +
geom_smooth(data=w2b_table, method="loess",
mapping = aes(x = dates, y = wing2body), colour="black") +
geom_point(data=w2b_table, mapping = aes(x = dates, y = wing2body)) +
ylim(0.71, 0.75) +
stat_smooth(data=last_n_years, mapping = aes(x = dates, y = wing2body), method = "lm", formula = y ~ x, se = FALSE,
colour="blue") # linetype = "longdash",
p0 = p0 + guides(fill = guide_legend(reverse = TRUE)) + theme_classic() +
theme(axis.text=element_text(size=13),
axis.title=element_text(size=16),
plot.title=element_text(size=20)) + guides(color = FALSE) + scale_linetype(guide = FALSE) +
theme(legend.key = element_rect(fill = "white", color = NA),
legend.key.size = unit(0.5, "cm"),
legend.key.width = unit(0.5,"cm")) +
scale_fill_manual(values = "blue", labels=c("C. corindum", "K. elegans")) +
labs(fill = "Host Plant") + theme(legend.title = element_text(size = 14),
legend.text = element_text(size = 13, face="italic")) +
scale_color_manual(values="blue") + theme(legend.position = c(0.2, 0.88)) + theme(legend.title = element_blank())
# loess
alpha = paste("alpha[loess]==", ggplot_build(p0)$data[[2]]$alpha[1])## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'degree="lambda[loess]==0"
mlinear = glm(wing2body ~ dates, data=d, family=gaussian)
pvalue = paste0("italic(p)[glm]==",round(summary(mlinear)$coeff[[8]],2))
p0 = p0 + annotate(geom="text", x=unique(d$dates)[10], y=0.74, label=alpha, color="black", parse=TRUE, size=6) +
annotate(geom="text", x=unique(d$dates)[10], y=0.748, label=degree, color="black",parse=TRUE, size=6) +
annotate(geom="text", x=unique(d$dates)[15], y=0.718, label=pvalue, color="black", parse=TRUE, size=6)
# last 5 years
mllinear = glm(wing2body ~ dates, data=w2b_table[w2b_table$date_b=="2015-2020",], family=gaussian)
pval = paste0("italic(p)[glm]==",round(summary(mllinear)$coeff[[8]],2))
summary(mllinear)##
## Call:
## glm(formula = wing2body ~ dates, family = gaussian, data = w2b_table[w2b_table$date_b ==
## "2015-2020", ])
##
## Deviance Residuals:
## 5 6 7 8 9
## -0.0009906 -0.0010888 0.0046361 0.0015432 -0.0040999
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.594e-01 9.356e-02 9.185 0.00273 **
## dates -7.209e-06 5.218e-06 -1.381 0.26106
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 1.428347e-05)
##
## Null deviance: 7.0107e-05 on 4 degrees of freedom
## Residual deviance: 4.2850e-05 on 3 degrees of freedom
## AIC: -38.147
##
## Number of Fisher Scoring iterations: 2p0 = p0 + annotate(geom="text", x=unique(d$dates)[24], y=0.743, label=pval, color="blue", parse=TRUE, size=6)
p0## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
head(ggplot_build(p0)$data)## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'## [[1]]
## xintercept PANEL group colour size linetype alpha
## 1 2013-04-01 1 -1 gainsboro 0.5 1 NA
## 2 2015-04-01 1 -1 gainsboro 0.5 1 NA
## 3 2016-12-01 1 -1 gainsboro 0.5 1 NA
## 4 2017-08-01 1 -1 gainsboro 0.5 1 NA
## 5 2018-09-01 1 -1 gainsboro 0.5 1 NA
## 6 2019-05-01 1 -1 gainsboro 0.5 1 NA
## 7 2019-10-01 1 -1 gainsboro 0.5 1 NA
## 8 2020-02-01 1 -1 gainsboro 0.5 1 NA
##
## [[2]]
## x y flipped_aes PANEL group colour fill size linetype
## 1 15796.00 0.7309072 FALSE 1 -1 black grey60 0.5 dashed
## 2 15827.61 0.7308962 FALSE 1 -1 black grey60 0.5 dashed
## 3 15859.22 0.7308853 FALSE 1 -1 black grey60 0.5 dashed
## 4 15890.82 0.7308744 FALSE 1 -1 black grey60 0.5 dashed
## 5 15922.43 0.7308635 FALSE 1 -1 black grey60 0.5 dashed
## 6 15954.04 0.7308525 FALSE 1 -1 black grey60 0.5 dashed
## 7 15985.65 0.7308416 FALSE 1 -1 black grey60 0.5 dashed
## 8 16017.25 0.7308307 FALSE 1 -1 black grey60 0.5 dashed
## 9 16048.86 0.7308198 FALSE 1 -1 black grey60 0.5 dashed
## 10 16080.47 0.7308088 FALSE 1 -1 black grey60 0.5 dashed
## 11 16112.08 0.7307979 FALSE 1 -1 black grey60 0.5 dashed
## 12 16143.68 0.7307870 FALSE 1 -1 black grey60 0.5 dashed
## 13 16175.29 0.7307761 FALSE 1 -1 black grey60 0.5 dashed
## 14 16206.90 0.7307651 FALSE 1 -1 black grey60 0.5 dashed
## 15 16238.51 0.7307542 FALSE 1 -1 black grey60 0.5 dashed
## 16 16270.11 0.7307433 FALSE 1 -1 black grey60 0.5 dashed
## 17 16301.72 0.7307324 FALSE 1 -1 black grey60 0.5 dashed
## 18 16333.33 0.7307214 FALSE 1 -1 black grey60 0.5 dashed
## 19 16364.94 0.7307105 FALSE 1 -1 black grey60 0.5 dashed
## 20 16396.54 0.7306996 FALSE 1 -1 black grey60 0.5 dashed
## 21 16428.15 0.7306887 FALSE 1 -1 black grey60 0.5 dashed
## 22 16459.76 0.7306777 FALSE 1 -1 black grey60 0.5 dashed
## 23 16491.37 0.7306668 FALSE 1 -1 black grey60 0.5 dashed
## 24 16522.97 0.7306559 FALSE 1 -1 black grey60 0.5 dashed
## 25 16554.58 0.7306450 FALSE 1 -1 black grey60 0.5 dashed
## 26 16586.19 0.7306340 FALSE 1 -1 black grey60 0.5 dashed
## 27 16617.80 0.7306231 FALSE 1 -1 black grey60 0.5 dashed
## 28 16649.41 0.7306122 FALSE 1 -1 black grey60 0.5 dashed
## 29 16681.01 0.7306013 FALSE 1 -1 black grey60 0.5 dashed
## 30 16712.62 0.7305903 FALSE 1 -1 black grey60 0.5 dashed
## 31 16744.23 0.7305794 FALSE 1 -1 black grey60 0.5 dashed
## 32 16775.84 0.7305685 FALSE 1 -1 black grey60 0.5 dashed
## 33 16807.44 0.7305576 FALSE 1 -1 black grey60 0.5 dashed
## 34 16839.05 0.7305466 FALSE 1 -1 black grey60 0.5 dashed
## 35 16870.66 0.7305357 FALSE 1 -1 black grey60 0.5 dashed
## 36 16902.27 0.7305248 FALSE 1 -1 black grey60 0.5 dashed
## 37 16933.87 0.7305139 FALSE 1 -1 black grey60 0.5 dashed
## 38 16965.48 0.7305029 FALSE 1 -1 black grey60 0.5 dashed
## 39 16997.09 0.7304920 FALSE 1 -1 black grey60 0.5 dashed
## 40 17028.70 0.7304811 FALSE 1 -1 black grey60 0.5 dashed
## 41 17060.30 0.7304702 FALSE 1 -1 black grey60 0.5 dashed
## 42 17091.91 0.7304592 FALSE 1 -1 black grey60 0.5 dashed
## 43 17123.52 0.7304483 FALSE 1 -1 black grey60 0.5 dashed
## 44 17155.13 0.7304374 FALSE 1 -1 black grey60 0.5 dashed
## 45 17186.73 0.7304265 FALSE 1 -1 black grey60 0.5 dashed
## 46 17218.34 0.7304155 FALSE 1 -1 black grey60 0.5 dashed
## 47 17249.95 0.7304046 FALSE 1 -1 black grey60 0.5 dashed
## 48 17281.56 0.7303937 FALSE 1 -1 black grey60 0.5 dashed
## 49 17313.16 0.7303828 FALSE 1 -1 black grey60 0.5 dashed
## 50 17344.77 0.7303718 FALSE 1 -1 black grey60 0.5 dashed
## 51 17376.38 0.7303609 FALSE 1 -1 black grey60 0.5 dashed
## 52 17407.99 0.7303500 FALSE 1 -1 black grey60 0.5 dashed
## 53 17439.59 0.7303391 FALSE 1 -1 black grey60 0.5 dashed
## 54 17471.20 0.7303281 FALSE 1 -1 black grey60 0.5 dashed
## 55 17502.81 0.7303172 FALSE 1 -1 black grey60 0.5 dashed
## 56 17534.42 0.7303063 FALSE 1 -1 black grey60 0.5 dashed
## 57 17566.03 0.7302954 FALSE 1 -1 black grey60 0.5 dashed
## 58 17597.63 0.7302844 FALSE 1 -1 black grey60 0.5 dashed
## 59 17629.24 0.7302735 FALSE 1 -1 black grey60 0.5 dashed
## 60 17660.85 0.7302626 FALSE 1 -1 black grey60 0.5 dashed
## 61 17692.46 0.7302517 FALSE 1 -1 black grey60 0.5 dashed
## 62 17724.06 0.7302407 FALSE 1 -1 black grey60 0.5 dashed
## 63 17755.67 0.7302298 FALSE 1 -1 black grey60 0.5 dashed
## 64 17787.28 0.7302189 FALSE 1 -1 black grey60 0.5 dashed
## 65 17818.89 0.7302080 FALSE 1 -1 black grey60 0.5 dashed
## 66 17850.49 0.7301970 FALSE 1 -1 black grey60 0.5 dashed
## 67 17882.10 0.7301861 FALSE 1 -1 black grey60 0.5 dashed
## 68 17913.71 0.7301752 FALSE 1 -1 black grey60 0.5 dashed
## 69 17945.32 0.7301643 FALSE 1 -1 black grey60 0.5 dashed
## 70 17976.92 0.7301533 FALSE 1 -1 black grey60 0.5 dashed
## 71 18008.53 0.7301424 FALSE 1 -1 black grey60 0.5 dashed
## 72 18040.14 0.7301315 FALSE 1 -1 black grey60 0.5 dashed
## 73 18071.75 0.7301206 FALSE 1 -1 black grey60 0.5 dashed
## 74 18103.35 0.7301096 FALSE 1 -1 black grey60 0.5 dashed
## 75 18134.96 0.7300987 FALSE 1 -1 black grey60 0.5 dashed
## 76 18166.57 0.7300878 FALSE 1 -1 black grey60 0.5 dashed
## 77 18198.18 0.7300769 FALSE 1 -1 black grey60 0.5 dashed
## 78 18229.78 0.7300659 FALSE 1 -1 black grey60 0.5 dashed
## 79 18261.39 0.7300550 FALSE 1 -1 black grey60 0.5 dashed
## 80 18293.00 0.7300441 FALSE 1 -1 black grey60 0.5 dashed
## weight alpha
## 1 1 0.4
## 2 1 0.4
## 3 1 0.4
## 4 1 0.4
## 5 1 0.4
## 6 1 0.4
## 7 1 0.4
## 8 1 0.4
## 9 1 0.4
## 10 1 0.4
## 11 1 0.4
## 12 1 0.4
## 13 1 0.4
## 14 1 0.4
## 15 1 0.4
## 16 1 0.4
## 17 1 0.4
## 18 1 0.4
## 19 1 0.4
## 20 1 0.4
## 21 1 0.4
## 22 1 0.4
## 23 1 0.4
## 24 1 0.4
## 25 1 0.4
## 26 1 0.4
## 27 1 0.4
## 28 1 0.4
## 29 1 0.4
## 30 1 0.4
## 31 1 0.4
## 32 1 0.4
## 33 1 0.4
## 34 1 0.4
## 35 1 0.4
## 36 1 0.4
## 37 1 0.4
## 38 1 0.4
## 39 1 0.4
## 40 1 0.4
## 41 1 0.4
## 42 1 0.4
## 43 1 0.4
## 44 1 0.4
## 45 1 0.4
## 46 1 0.4
## 47 1 0.4
## 48 1 0.4
## 49 1 0.4
## 50 1 0.4
## 51 1 0.4
## 52 1 0.4
## 53 1 0.4
## 54 1 0.4
## 55 1 0.4
## 56 1 0.4
## 57 1 0.4
## 58 1 0.4
## 59 1 0.4
## 60 1 0.4
## 61 1 0.4
## 62 1 0.4
## 63 1 0.4
## 64 1 0.4
## 65 1 0.4
## 66 1 0.4
## 67 1 0.4
## 68 1 0.4
## 69 1 0.4
## 70 1 0.4
## 71 1 0.4
## 72 1 0.4
## 73 1 0.4
## 74 1 0.4
## 75 1 0.4
## 76 1 0.4
## 77 1 0.4
## 78 1 0.4
## 79 1 0.4
## 80 1 0.4
##
## [[3]]
## x y ymin ymax se flipped_aes PANEL group
## 1 15796.00 0.7321013 0.7221142 0.7420883 0.003243019 FALSE 1 -1
## 2 15827.61 0.7316196 0.7223835 0.7408557 0.002999189 FALSE 1 -1
## 3 15859.22 0.7311616 0.7225927 0.7397305 0.002782527 FALSE 1 -1
## 4 15890.82 0.7307275 0.7227373 0.7387176 0.002594589 FALSE 1 -1
## 5 15922.43 0.7303172 0.7228132 0.7378212 0.002436730 FALSE 1 -1
## 6 15954.04 0.7299310 0.7228177 0.7370442 0.002309846 FALSE 1 -1
## 7 15985.65 0.7295689 0.7227506 0.7363872 0.002214074 FALSE 1 -1
## 8 16017.25 0.7292310 0.7226145 0.7358475 0.002148540 FALSE 1 -1
## 9 16048.86 0.7289175 0.7224158 0.7354192 0.002111249 FALSE 1 -1
## 10 16080.47 0.7286285 0.7221640 0.7350929 0.002099171 FALSE 1 -1
## 11 16112.08 0.7283640 0.7218707 0.7348572 0.002108521 FALSE 1 -1
## 12 16143.68 0.7281241 0.7215489 0.7346994 0.002135137 FALSE 1 -1
## 13 16175.29 0.7279096 0.7212127 0.7346065 0.002174647 FALSE 1 -1
## 14 16206.90 0.7277238 0.7208820 0.7345657 0.002221709 FALSE 1 -1
## 15 16238.51 0.7275665 0.7205672 0.7345658 0.002272842 FALSE 1 -1
## 16 16270.11 0.7274364 0.7202747 0.7345981 0.002325582 FALSE 1 -1
## 17 16301.72 0.7273325 0.7200093 0.7346558 0.002378043 FALSE 1 -1
## 18 16333.33 0.7272537 0.7197742 0.7347332 0.002428773 FALSE 1 -1
## 19 16364.94 0.7271987 0.7195719 0.7348256 0.002476632 FALSE 1 -1
## 20 16396.54 0.7271665 0.7194039 0.7349291 0.002520703 FALSE 1 -1
## 21 16428.15 0.7271560 0.7192717 0.7350403 0.002560218 FALSE 1 -1
## 22 16459.76 0.7271659 0.7191760 0.7351558 0.002594516 FALSE 1 -1
## 23 16491.37 0.7271952 0.7191175 0.7352729 0.002623019 FALSE 1 -1
## 24 16522.97 0.7272427 0.7190966 0.7353887 0.002645221 FALSE 1 -1
## 25 16554.58 0.7273393 0.7191696 0.7355089 0.002652885 FALSE 1 -1
## 26 16586.19 0.7275246 0.7194037 0.7356454 0.002637029 FALSE 1 -1
## 27 16617.80 0.7277892 0.7197759 0.7358026 0.002602115 FALSE 1 -1
## 28 16649.41 0.7281236 0.7202603 0.7359869 0.002553395 FALSE 1 -1
## 29 16681.01 0.7285179 0.7208295 0.7362063 0.002496607 FALSE 1 -1
## 30 16712.62 0.7289626 0.7214554 0.7364697 0.002437750 FALSE 1 -1
## 31 16744.23 0.7294478 0.7221100 0.7367857 0.002382783 FALSE 1 -1
## 32 16775.84 0.7299640 0.7227664 0.7371615 0.002337214 FALSE 1 -1
## 33 16807.44 0.7305013 0.7234011 0.7376015 0.002305604 FALSE 1 -1
## 34 16839.05 0.7310502 0.7239949 0.7381055 0.002291028 FALSE 1 -1
## 35 16870.66 0.7316009 0.7245344 0.7386674 0.002294660 FALSE 1 -1
## 36 16902.27 0.7321437 0.7250128 0.7392747 0.002315589 FALSE 1 -1
## 37 16933.87 0.7326690 0.7254291 0.7399088 0.002350952 FALSE 1 -1
## 38 16965.48 0.7331670 0.7257874 0.7405466 0.002396328 FALSE 1 -1
## 39 16997.09 0.7336280 0.7260946 0.7411614 0.002446274 FALSE 1 -1
## 40 17028.70 0.7340424 0.7263594 0.7417254 0.002494850 FALSE 1 -1
## 41 17060.30 0.7344004 0.7265905 0.7422104 0.002536068 FALSE 1 -1
## 42 17091.91 0.7346924 0.7267958 0.7425891 0.002564229 FALSE 1 -1
## 43 17123.52 0.7349087 0.7269814 0.7428361 0.002574195 FALSE 1 -1
## 44 17155.13 0.7350202 0.7271411 0.7428993 0.002558532 FALSE 1 -1
## 45 17186.73 0.7349584 0.7272259 0.7426909 0.002510930 FALSE 1 -1
## 46 17218.34 0.7347490 0.7271967 0.7423012 0.002452399 FALSE 1 -1
## 47 17249.95 0.7344272 0.7270102 0.7418442 0.002408482 FALSE 1 -1
## 48 17281.56 0.7340282 0.7266337 0.7414227 0.002401162 FALSE 1 -1
## 49 17313.16 0.7335872 0.7260666 0.7411078 0.002442126 FALSE 1 -1
## 50 17344.77 0.7331393 0.7253523 0.7409263 0.002528635 FALSE 1 -1
## 51 17376.38 0.7327197 0.7245753 0.7408641 0.002644675 FALSE 1 -1
## 52 17407.99 0.7323773 0.7239211 0.7408334 0.002745907 FALSE 1 -1
## 53 17439.59 0.7321357 0.7235313 0.7407401 0.002794057 FALSE 1 -1
## 54 17471.20 0.7319811 0.7233695 0.7405928 0.002796400 FALSE 1 -1
## 55 17502.81 0.7318993 0.7233775 0.7404212 0.002767245 FALSE 1 -1
## 56 17534.42 0.7318762 0.7234879 0.7402645 0.002723874 FALSE 1 -1
## 57 17566.03 0.7318976 0.7236337 0.7401614 0.002683467 FALSE 1 -1
## 58 17597.63 0.7319494 0.7237579 0.7401409 0.002659963 FALSE 1 -1
## 59 17629.24 0.7320175 0.7238227 0.7402123 0.002661056 FALSE 1 -1
## 60 17660.85 0.7320878 0.7238153 0.7403602 0.002686271 FALSE 1 -1
## 61 17692.46 0.7321460 0.7237481 0.7405440 0.002727022 FALSE 1 -1
## 62 17724.06 0.7321782 0.7236522 0.7407043 0.002768604 FALSE 1 -1
## 63 17755.67 0.7321702 0.7235686 0.7407718 0.002793143 FALSE 1 -1
## 64 17787.28 0.7321183 0.7235487 0.7406880 0.002782755 FALSE 1 -1
## 65 17818.89 0.7321054 0.7236920 0.7405188 0.002732026 FALSE 1 -1
## 66 17850.49 0.7321246 0.7239545 0.7402947 0.002653025 FALSE 1 -1
## 67 17882.10 0.7321418 0.7242620 0.7400216 0.002558763 FALSE 1 -1
## 68 17913.71 0.7321231 0.7245502 0.7396959 0.002459076 FALSE 1 -1
## 69 17945.32 0.7320341 0.7247703 0.7392980 0.002358748 FALSE 1 -1
## 70 17976.92 0.7318410 0.7248896 0.7387924 0.002257293 FALSE 1 -1
## 71 18008.53 0.7315095 0.7248844 0.7381345 0.002151314 FALSE 1 -1
## 72 18040.14 0.7310496 0.7247592 0.7373399 0.002042632 FALSE 1 -1
## 73 18071.75 0.7305501 0.7245566 0.7365435 0.001946226 FALSE 1 -1
## 74 18103.35 0.7299995 0.7242167 0.7357823 0.001877812 FALSE 1 -1
## 75 18134.96 0.7293808 0.7236648 0.7350968 0.001856127 FALSE 1 -1
## 76 18166.57 0.7286767 0.7228235 0.7345299 0.001900676 FALSE 1 -1
## 77 18198.18 0.7278815 0.7216367 0.7341263 0.002027837 FALSE 1 -1
## 78 18229.78 0.7270075 0.7200921 0.7339230 0.002245614 FALSE 1 -1
## 79 18261.39 0.7260586 0.7182007 0.7339165 0.002551653 FALSE 1 -1
## 80 18293.00 0.7250384 0.7159909 0.7340859 0.002937927 FALSE 1 -1
## colour fill size linetype weight alpha
## 1 black grey60 1 1 1 0.4
## 2 black grey60 1 1 1 0.4
## 3 black grey60 1 1 1 0.4
## 4 black grey60 1 1 1 0.4
## 5 black grey60 1 1 1 0.4
## 6 black grey60 1 1 1 0.4
## 7 black grey60 1 1 1 0.4
## 8 black grey60 1 1 1 0.4
## 9 black grey60 1 1 1 0.4
## 10 black grey60 1 1 1 0.4
## 11 black grey60 1 1 1 0.4
## 12 black grey60 1 1 1 0.4
## 13 black grey60 1 1 1 0.4
## 14 black grey60 1 1 1 0.4
## 15 black grey60 1 1 1 0.4
## 16 black grey60 1 1 1 0.4
## 17 black grey60 1 1 1 0.4
## 18 black grey60 1 1 1 0.4
## 19 black grey60 1 1 1 0.4
## 20 black grey60 1 1 1 0.4
## 21 black grey60 1 1 1 0.4
## 22 black grey60 1 1 1 0.4
## 23 black grey60 1 1 1 0.4
## 24 black grey60 1 1 1 0.4
## 25 black grey60 1 1 1 0.4
## 26 black grey60 1 1 1 0.4
## 27 black grey60 1 1 1 0.4
## 28 black grey60 1 1 1 0.4
## 29 black grey60 1 1 1 0.4
## 30 black grey60 1 1 1 0.4
## 31 black grey60 1 1 1 0.4
## 32 black grey60 1 1 1 0.4
## 33 black grey60 1 1 1 0.4
## 34 black grey60 1 1 1 0.4
## 35 black grey60 1 1 1 0.4
## 36 black grey60 1 1 1 0.4
## 37 black grey60 1 1 1 0.4
## 38 black grey60 1 1 1 0.4
## 39 black grey60 1 1 1 0.4
## 40 black grey60 1 1 1 0.4
## 41 black grey60 1 1 1 0.4
## 42 black grey60 1 1 1 0.4
## 43 black grey60 1 1 1 0.4
## 44 black grey60 1 1 1 0.4
## 45 black grey60 1 1 1 0.4
## 46 black grey60 1 1 1 0.4
## 47 black grey60 1 1 1 0.4
## 48 black grey60 1 1 1 0.4
## 49 black grey60 1 1 1 0.4
## 50 black grey60 1 1 1 0.4
## 51 black grey60 1 1 1 0.4
## 52 black grey60 1 1 1 0.4
## 53 black grey60 1 1 1 0.4
## 54 black grey60 1 1 1 0.4
## 55 black grey60 1 1 1 0.4
## 56 black grey60 1 1 1 0.4
## 57 black grey60 1 1 1 0.4
## 58 black grey60 1 1 1 0.4
## 59 black grey60 1 1 1 0.4
## 60 black grey60 1 1 1 0.4
## 61 black grey60 1 1 1 0.4
## 62 black grey60 1 1 1 0.4
## 63 black grey60 1 1 1 0.4
## 64 black grey60 1 1 1 0.4
## 65 black grey60 1 1 1 0.4
## 66 black grey60 1 1 1 0.4
## 67 black grey60 1 1 1 0.4
## 68 black grey60 1 1 1 0.4
## 69 black grey60 1 1 1 0.4
## 70 black grey60 1 1 1 0.4
## 71 black grey60 1 1 1 0.4
## 72 black grey60 1 1 1 0.4
## 73 black grey60 1 1 1 0.4
## 74 black grey60 1 1 1 0.4
## 75 black grey60 1 1 1 0.4
## 76 black grey60 1 1 1 0.4
## 77 black grey60 1 1 1 0.4
## 78 black grey60 1 1 1 0.4
## 79 black grey60 1 1 1 0.4
## 80 black grey60 1 1 1 0.4
##
## [[4]]
## y x PANEL group shape colour size fill alpha stroke
## 1 0.7313394 15796 1 -1 19 black 1.5 NA NA 0.5
## 2 0.7304482 16161 1 -1 19 black 1.5 NA NA 0.5
## 3 0.7243691 16526 1 -1 19 black 1.5 NA NA 0.5
## 4 0.7366048 17136 1 -1 19 black 1.5 NA NA 0.5
## 5 0.7331342 17379 1 -1 19 black 1.5 NA NA 0.5
## 6 0.7301813 17775 1 -1 19 black 1.5 NA NA 0.5
## 7 0.7341616 18017 1 -1 19 black 1.5 NA NA 0.5
## 8 0.7299658 18170 1 -1 19 black 1.5 NA NA 0.5
## 9 0.7234360 18293 1 -1 19 black 1.5 NA NA 0.5
##
## [[5]]
## x y flipped_aes PANEL group colour fill size linetype
## 1 17379.00 0.7341247 FALSE 1 -1 blue grey60 1 1
## 2 17390.57 0.7340413 FALSE 1 -1 blue grey60 1 1
## 3 17402.14 0.7339579 FALSE 1 -1 blue grey60 1 1
## 4 17413.71 0.7338745 FALSE 1 -1 blue grey60 1 1
## 5 17425.28 0.7337911 FALSE 1 -1 blue grey60 1 1
## 6 17436.85 0.7337077 FALSE 1 -1 blue grey60 1 1
## 7 17448.42 0.7336243 FALSE 1 -1 blue grey60 1 1
## 8 17459.99 0.7335409 FALSE 1 -1 blue grey60 1 1
## 9 17471.56 0.7334575 FALSE 1 -1 blue grey60 1 1
## 10 17483.13 0.7333741 FALSE 1 -1 blue grey60 1 1
## 11 17494.70 0.7332907 FALSE 1 -1 blue grey60 1 1
## 12 17506.27 0.7332073 FALSE 1 -1 blue grey60 1 1
## 13 17517.84 0.7331239 FALSE 1 -1 blue grey60 1 1
## 14 17529.41 0.7330405 FALSE 1 -1 blue grey60 1 1
## 15 17540.97 0.7329571 FALSE 1 -1 blue grey60 1 1
## 16 17552.54 0.7328737 FALSE 1 -1 blue grey60 1 1
## 17 17564.11 0.7327903 FALSE 1 -1 blue grey60 1 1
## 18 17575.68 0.7327069 FALSE 1 -1 blue grey60 1 1
## 19 17587.25 0.7326235 FALSE 1 -1 blue grey60 1 1
## 20 17598.82 0.7325401 FALSE 1 -1 blue grey60 1 1
## 21 17610.39 0.7324567 FALSE 1 -1 blue grey60 1 1
## 22 17621.96 0.7323733 FALSE 1 -1 blue grey60 1 1
## 23 17633.53 0.7322899 FALSE 1 -1 blue grey60 1 1
## 24 17645.10 0.7322065 FALSE 1 -1 blue grey60 1 1
## 25 17656.67 0.7321231 FALSE 1 -1 blue grey60 1 1
## 26 17668.24 0.7320397 FALSE 1 -1 blue grey60 1 1
## 27 17679.81 0.7319563 FALSE 1 -1 blue grey60 1 1
## 28 17691.38 0.7318729 FALSE 1 -1 blue grey60 1 1
## 29 17702.95 0.7317895 FALSE 1 -1 blue grey60 1 1
## 30 17714.52 0.7317060 FALSE 1 -1 blue grey60 1 1
## 31 17726.09 0.7316226 FALSE 1 -1 blue grey60 1 1
## 32 17737.66 0.7315392 FALSE 1 -1 blue grey60 1 1
## 33 17749.23 0.7314558 FALSE 1 -1 blue grey60 1 1
## 34 17760.80 0.7313724 FALSE 1 -1 blue grey60 1 1
## 35 17772.37 0.7312890 FALSE 1 -1 blue grey60 1 1
## 36 17783.94 0.7312056 FALSE 1 -1 blue grey60 1 1
## 37 17795.51 0.7311222 FALSE 1 -1 blue grey60 1 1
## 38 17807.08 0.7310388 FALSE 1 -1 blue grey60 1 1
## 39 17818.65 0.7309554 FALSE 1 -1 blue grey60 1 1
## 40 17830.22 0.7308720 FALSE 1 -1 blue grey60 1 1
## 41 17841.78 0.7307886 FALSE 1 -1 blue grey60 1 1
## 42 17853.35 0.7307052 FALSE 1 -1 blue grey60 1 1
## 43 17864.92 0.7306218 FALSE 1 -1 blue grey60 1 1
## 44 17876.49 0.7305384 FALSE 1 -1 blue grey60 1 1
## 45 17888.06 0.7304550 FALSE 1 -1 blue grey60 1 1
## 46 17899.63 0.7303716 FALSE 1 -1 blue grey60 1 1
## 47 17911.20 0.7302882 FALSE 1 -1 blue grey60 1 1
## 48 17922.77 0.7302048 FALSE 1 -1 blue grey60 1 1
## 49 17934.34 0.7301214 FALSE 1 -1 blue grey60 1 1
## 50 17945.91 0.7300380 FALSE 1 -1 blue grey60 1 1
## 51 17957.48 0.7299546 FALSE 1 -1 blue grey60 1 1
## 52 17969.05 0.7298712 FALSE 1 -1 blue grey60 1 1
## 53 17980.62 0.7297878 FALSE 1 -1 blue grey60 1 1
## 54 17992.19 0.7297044 FALSE 1 -1 blue grey60 1 1
## 55 18003.76 0.7296210 FALSE 1 -1 blue grey60 1 1
## 56 18015.33 0.7295376 FALSE 1 -1 blue grey60 1 1
## 57 18026.90 0.7294542 FALSE 1 -1 blue grey60 1 1
## 58 18038.47 0.7293708 FALSE 1 -1 blue grey60 1 1
## 59 18050.04 0.7292874 FALSE 1 -1 blue grey60 1 1
## 60 18061.61 0.7292040 FALSE 1 -1 blue grey60 1 1
## 61 18073.18 0.7291206 FALSE 1 -1 blue grey60 1 1
## 62 18084.75 0.7290372 FALSE 1 -1 blue grey60 1 1
## 63 18096.32 0.7289538 FALSE 1 -1 blue grey60 1 1
## 64 18107.89 0.7288704 FALSE 1 -1 blue grey60 1 1
## 65 18119.46 0.7287870 FALSE 1 -1 blue grey60 1 1
## 66 18131.03 0.7287036 FALSE 1 -1 blue grey60 1 1
## 67 18142.59 0.7286202 FALSE 1 -1 blue grey60 1 1
## 68 18154.16 0.7285368 FALSE 1 -1 blue grey60 1 1
## 69 18165.73 0.7284534 FALSE 1 -1 blue grey60 1 1
## 70 18177.30 0.7283700 FALSE 1 -1 blue grey60 1 1
## 71 18188.87 0.7282866 FALSE 1 -1 blue grey60 1 1
## 72 18200.44 0.7282032 FALSE 1 -1 blue grey60 1 1
## 73 18212.01 0.7281198 FALSE 1 -1 blue grey60 1 1
## 74 18223.58 0.7280364 FALSE 1 -1 blue grey60 1 1
## 75 18235.15 0.7279530 FALSE 1 -1 blue grey60 1 1
## 76 18246.72 0.7278695 FALSE 1 -1 blue grey60 1 1
## 77 18258.29 0.7277861 FALSE 1 -1 blue grey60 1 1
## 78 18269.86 0.7277027 FALSE 1 -1 blue grey60 1 1
## 79 18281.43 0.7276193 FALSE 1 -1 blue grey60 1 1
## 80 18293.00 0.7275359 FALSE 1 -1 blue grey60 1 1
## weight alpha
## 1 1 0.4
## 2 1 0.4
## 3 1 0.4
## 4 1 0.4
## 5 1 0.4
## 6 1 0.4
## 7 1 0.4
## 8 1 0.4
## 9 1 0.4
## 10 1 0.4
## 11 1 0.4
## 12 1 0.4
## 13 1 0.4
## 14 1 0.4
## 15 1 0.4
## 16 1 0.4
## 17 1 0.4
## 18 1 0.4
## 19 1 0.4
## 20 1 0.4
## 21 1 0.4
## 22 1 0.4
## 23 1 0.4
## 24 1 0.4
## 25 1 0.4
## 26 1 0.4
## 27 1 0.4
## 28 1 0.4
## 29 1 0.4
## 30 1 0.4
## 31 1 0.4
## 32 1 0.4
## 33 1 0.4
## 34 1 0.4
## 35 1 0.4
## 36 1 0.4
## 37 1 0.4
## 38 1 0.4
## 39 1 0.4
## 40 1 0.4
## 41 1 0.4
## 42 1 0.4
## 43 1 0.4
## 44 1 0.4
## 45 1 0.4
## 46 1 0.4
## 47 1 0.4
## 48 1 0.4
## 49 1 0.4
## 50 1 0.4
## 51 1 0.4
## 52 1 0.4
## 53 1 0.4
## 54 1 0.4
## 55 1 0.4
## 56 1 0.4
## 57 1 0.4
## 58 1 0.4
## 59 1 0.4
## 60 1 0.4
## 61 1 0.4
## 62 1 0.4
## 63 1 0.4
## 64 1 0.4
## 65 1 0.4
## 66 1 0.4
## 67 1 0.4
## 68 1 0.4
## 69 1 0.4
## 70 1 0.4
## 71 1 0.4
## 72 1 0.4
## 73 1 0.4
## 74 1 0.4
## 75 1 0.4
## 76 1 0.4
## 77 1 0.4
## 78 1 0.4
## 79 1 0.4
## 80 1 0.4
##
## [[6]]
## y x label PANEL group colour size angle hjust vjust
## 1 0.74 16526.05 alpha[loess]== 0.4 1 -1 black 6 0 0.5 0.5
## alpha family fontface lineheight
## 1 NA 1 1.2w2b_summary<-aggregate(wing2body~sex*pophost*month_of_year, data=data_long, FUN=mean)
df = w2b_summary
xlab_dates = na.omit(sort(unique(data_long$month_of_year)))
xlab_months = xlab_dates[c(-2,-5)]
month_labs <- c("Feb", "May", "Aug", "Oct", "Dec")
sex_colors_shade = c("brown1", "sienna4")
sex_colors_pts = c("brown1", "grey27")df$pophost[df$pophost=="C.corindum"]<-"C. corindum"
df$pophost[df$pophost=="K.elegans"]<-"K. elegans"
df$pophost = factor(df$pophost, levels = c("K. elegans", "C. corindum") )
df$`Host Plant` = df$pophost
df$sex[df$sex=="F"]<-"Females"
df$sex[df$sex=="M"]<-"Males"
df$sex = factor(df$sex, levels = c("Males", "Females") )
df$`Sex` = df$sex customPlot = list( theme_classic(),
theme(axis.text=element_text(size=13),
axis.title=element_text(size=16),
plot.title=element_text(size=20),),
theme(legend.position = c(0.2, 0.9)),
theme(legend.title = element_text(size=14, face="italic"),
legend.text = element_text(size = 13, face="italic"))
)# Females
p1 = ggplot() +
ggtitle("Females") + xlab("Year") + ylab("Wing-to-Body Ratio") + # title and labels
geom_vline(xintercept = xlab_years, color="gainsboro") + # add grey lines for the each collection date
geom_smooth(data=dfF, method="glm", se=FALSE, linetype = "dashed", # GLM regression and pvalue
mapping = aes(x = dates, y = wing2body), colour="black", lwd=0.5) +
geom_smooth(data=dfF, method="loess", # LOESS regression
mapping = aes(x = dates, y = wing2body, colour=`Host Plant`, fill=`Host Plant`)) +
geom_point(data=dfF, mapping = aes(x = dates, y = wing2body, colour=`Host Plant`)) + # data points and yaxis value range
ylim(0.70, 0.76) +
customPlot +
scale_color_manual(values=c("C. corindum" = "turquoise3", "K. elegans" = "springgreen4")) +
scale_fill_manual(values = c("C. corindum" = "turquoise3", "K. elegans" = "green"))
# Males
p2 = ggplot() +
ggtitle("Males") + xlab("Year") + ylab("Wing-to-Body Ratio") + # title and labels
geom_vline(xintercept = xlab_years, color="gainsboro") + # add grey lines for the each collection date
geom_smooth(data=dfM, method="glm", se=FALSE, linetype = "dashed", # GLM regression and pvalue
mapping = aes(x = dates, y = wing2body), colour="black", lwd=0.5) +
geom_smooth(data=dfM, method="loess", # LOESS regression
mapping = aes(x = dates, y = wing2body, colour=`Host Plant`, fill=`Host Plant`)) +
geom_point(data=dfM, mapping = aes(x = dates, y = wing2body, colour=`Host Plant`)) + # data points and yaxis value range
ylim(0.70, 0.76) +
customPlot +
scale_color_manual(values=c("C. corindum" = "turquoise3", "K. elegans" = "springgreen4")) +
scale_fill_manual(values = c("C. corindum" = "turquoise3", "K. elegans" = "green"))
grid.arrange(p1, p2, ncol=2)## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
Now do it by season because that’s more interesting…
p5 = ggplot() +
ggtitle("A") + xlab("Month") + ylab("Wing-to-Body Ratio") + # title and labels
geom_vline(xintercept = xlab_dates, color="gainsboro") + # add grey lines for the each collection date
geom_smooth(data=df, method="glm", se=FALSE, linetype = "dashed", # GLM regression and pvalue
mapping = aes(x = month_of_year, y = wing2body), colour="black", lwd=0.5) +
geom_smooth(data=df, method="loess", # LOESS regression
mapping = aes(x = month_of_year, y = wing2body, colour=`Host Plant`, fill=`Host Plant`)) +
geom_point(data=df, mapping = aes(x = month_of_year, y = wing2body, colour=`Host Plant`)) + # data points and yaxis value range
ylim(0.71, 0.75) +
customPlot +
scale_color_manual(values=c("C. corindum" = "turquoise3", "K. elegans" = "springgreen4")) +
scale_fill_manual(values = c("C. corindum" = "turquoise3", "K. elegans" = "green")) +
scale_x_continuous(breaks=xlab_months, labels= month_labs)
p6 = ggplot() +
ggtitle("B") + xlab("Month") + ylab("Wing-to-Body Ratio") + # title and labels
geom_vline(xintercept = xlab_dates, color="gainsboro") + # add grey lines for the each collection date
geom_smooth(data=df, method="glm", se=FALSE, linetype = "dashed", # GLM regression and pvalue
mapping = aes(x = month_of_year, y = wing2body), colour="black", lwd=0.5) +
geom_smooth(data=df, method="loess", # LOESS regression
mapping = aes(x = month_of_year, y = wing2body, colour=Sex, fill=Sex)) +
geom_point(data=df, mapping = aes(x = month_of_year, y = wing2body, colour=Sex)) + # data points and yaxis value range
ylim(0.71, 0.75) +
customPlot +
scale_color_manual(values=c("Females" = "brown3", "Males" = "black")) +
scale_fill_manual(values = c("Females" = "brown1", "Males" = "sienna4")) +
scale_x_continuous(breaks=xlab_months, labels= month_labs) + theme(axis.title.y = element_blank())
# extract the LOESS
alpha = paste("alpha[loess]==", ggplot_build(p5)$data[[2]]$alpha[1])## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'degree="lambda[loess]==0"
mlinear = glm(wing2body ~ month_of_year, data=df, family=gaussian)
#round(summary(mlinear)$coeff[[8]],3)
pvalue = paste0("italic(p)[glm]==",round(summary(mlinear)$coeff[[8]],3))
p5 = p5 + annotate(geom="text", x=10, y=0.744, label=alpha, color="black", parse=TRUE) +
annotate(geom="text", x=10, y=0.747, label=degree, color="black", parse=TRUE) +
annotate(geom="text", x=6.3, y=0.715, label=pvalue, color="black", parse=TRUE)
p6 = p6 + annotate(geom="text", x=10, y=0.744, label=alpha, color="black", parse=TRUE) +
annotate(geom="text", x=10, y=0.747, label=degree, color="black", parse=TRUE) +
annotate(geom="text", x=6.3, y=0.715, label=pvalue, color="black", parse=TRUE) figure = ggdraw() +
draw_plot(p5, 0, .5, .5, .5) +
draw_plot(p6, 0.5, .5, .5, .5) +
draw_plot(p0, 0, 0, 1, .5)## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'figure
#ggsave("w2b_over_time.pdf", plot=figure)Loess: a nonparametric, graphical tool for depicting relationships between variables. http://www.jtrive.com/docs/Jacoby_LOESS.pdf
# extract the LOESS
head(ggplot_build(p5)$data[[2]]) # alpha = 0.4 (smoothing parameter)## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'## x y flipped_aes PANEL group colour fill size linetype weight
## 1 2.000000 0.7252934 FALSE 1 -1 black grey60 0.5 dashed 1
## 2 2.126582 0.7253974 FALSE 1 -1 black grey60 0.5 dashed 1
## 3 2.253165 0.7255014 FALSE 1 -1 black grey60 0.5 dashed 1
## 4 2.379747 0.7256053 FALSE 1 -1 black grey60 0.5 dashed 1
## 5 2.506329 0.7257093 FALSE 1 -1 black grey60 0.5 dashed 1
## 6 2.632911 0.7258133 FALSE 1 -1 black grey60 0.5 dashed 1
## alpha
## 1 0.4
## 2 0.4
## 3 0.4
## 4 0.4
## 5 0.4
## 6 0.4head(ggplot_build(p6)$data[[2]]) # alpha = 0.4 (smoothing parameter); lambda = degree of the local polynomial ## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'## x y flipped_aes PANEL group colour fill size linetype weight
## 1 2.000000 0.7252934 FALSE 1 -1 black grey60 0.5 dashed 1
## 2 2.126582 0.7253974 FALSE 1 -1 black grey60 0.5 dashed 1
## 3 2.253165 0.7255014 FALSE 1 -1 black grey60 0.5 dashed 1
## 4 2.379747 0.7256053 FALSE 1 -1 black grey60 0.5 dashed 1
## 5 2.506329 0.7257093 FALSE 1 -1 black grey60 0.5 dashed 1
## 6 2.632911 0.7258133 FALSE 1 -1 black grey60 0.5 dashed 1
## alpha
## 1 0.4
## 2 0.4
## 3 0.4
## 4 0.4
## 5 0.4
## 6 0.4# degree is 0 (estimated below)
# Using a zero degree polynomial turns LOESS into a weighted moving average.
# degrees: https://www.itl.nist.gov/div898/handbook/pmd/section1/pmd144.htmTry ggdraw() http://www.sthda.com/english/wiki/wiki.php?id_contents=7930 Map separate lines: https://stats.idre.ucla.edu/r/faq/how-can-i-explore-different-smooths-in-ggplot2/
dfF = df[df$sex=="Females",]
dfM = df[df$sex=="Males",]
dfF$pophost = factor(dfF$pophost, levels = c("K. elegans", "C. corindum") )
dfF$`Host Plant` = dfF$pophost
dfM$pophost = factor(dfM$pophost, levels = c("K. elegans", "C. corindum") )
dfM$`Host Plant` = dfM$pophost# females
p3 = ggplot() +
ggtitle("Females") + xlab("Month") + ylab("Wing-to-Body Ratio") +
geom_vline(xintercept = xlab_dates, color="gainsboro") +
geom_smooth(data=dfF, method="lm", se=FALSE, linetype = "dashed", lwd=0.5,
mapping = aes(x = month_of_year, y = wing2body, colour = `Host Plant`)) +
geom_smooth(data=dfF, method="loess",
mapping = aes(x = month_of_year, y = wing2body, colour=`Host Plant`, fill=`Host Plant`)) +
geom_point(data=dfF, mapping = aes(x = month_of_year, y = wing2body, colour=`Host Plant`)) +
ylim(0.70, 0.765) +
customPlot +
scale_color_manual(values=c("C. corindum" = "turquoise3", "K. elegans" = "springgreen4")) +
scale_fill_manual(values = c("C. corindum" = "turquoise3", "K. elegans" = "green")) +
scale_x_continuous(breaks=xlab_months, labels= month_labs)
# males
p4 = ggplot() +
ggtitle("Males") + xlab("Month") + ylab("Wing-to-Body Ratio") +
geom_vline(xintercept = xlab_dates, color="gainsboro") +
geom_smooth(data=dfM, method="lm", se=FALSE, linetype = "dashed", lwd=0.5,
mapping = aes(x = month_of_year, y = wing2body, colour = `Host Plant`)) +
geom_smooth(data=dfM, method="loess",
mapping = aes(x = month_of_year, y = wing2body, colour=`Host Plant`, fill=`Host Plant`)) +
geom_point(data=dfM, mapping = aes(x = month_of_year, y = wing2body, colour=`Host Plant`)) +
ylim(0.70, 0.765) +
customPlot +
scale_color_manual(values=c("C. corindum" = "turquoise3", "K. elegans" = "springgreen4")) +
scale_fill_manual(values = c("C. corindum" = "turquoise3", "K. elegans" = "green")) +
scale_x_continuous(breaks=xlab_months, labels= month_labs)
grid.arrange(p3, p4, ncol=2)## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
# extract the LOESS
head(ggplot_build(p3)$data[[2]])## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'## colour x y flipped_aes PANEL group fill size linetype
## 1 springgreen4 2.000000 0.7289583 FALSE 1 1 grey60 0.5 dashed
## 2 springgreen4 2.126582 0.7290421 FALSE 1 1 grey60 0.5 dashed
## 3 springgreen4 2.253165 0.7291258 FALSE 1 1 grey60 0.5 dashed
## 4 springgreen4 2.379747 0.7292096 FALSE 1 1 grey60 0.5 dashed
## 5 springgreen4 2.506329 0.7292934 FALSE 1 1 grey60 0.5 dashed
## 6 springgreen4 2.632911 0.7293771 FALSE 1 1 grey60 0.5 dashed
## weight alpha
## 1 1 0.4
## 2 1 0.4
## 3 1 0.4
## 4 1 0.4
## 5 1 0.4
## 6 1 0.4head(ggplot_build(p4)$data[[2]])## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'## colour x y flipped_aes PANEL group fill size linetype
## 1 springgreen4 2.000000 0.7286637 FALSE 1 1 grey60 0.5 dashed
## 2 springgreen4 2.126582 0.7288033 FALSE 1 1 grey60 0.5 dashed
## 3 springgreen4 2.253165 0.7289429 FALSE 1 1 grey60 0.5 dashed
## 4 springgreen4 2.379747 0.7290826 FALSE 1 1 grey60 0.5 dashed
## 5 springgreen4 2.506329 0.7292222 FALSE 1 1 grey60 0.5 dashed
## 6 springgreen4 2.632911 0.7293618 FALSE 1 1 grey60 0.5 dashed
## weight alpha
## 1 1 0.4
## 2 1 0.4
## 3 1 0.4
## 4 1 0.4
## 5 1 0.4
## 6 1 0.4# these calculations give you similar tables. Here is an example:
mloess <- loess(wing2body ~ month_of_year, data=dfF, degree=0)
# Using a zero degree polynomial turns LOESS into a weighted moving average.
xrange <- range(dfF$month_of_year)
xseq <- seq(from=xrange[1], to=xrange[2], length=80)
pred <- predict(mloess, newdata = data.frame(month_of_year = xseq), se=TRUE)## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : degree must be at least
## 1 for vertex influence matrix
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : degree must be at least
## 1 for vertex influence matrix
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : degree must be at least
## 1 for vertex influence matrix
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : degree must be at least
## 1 for vertex influence matrix
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : degree must be at least
## 1 for vertex influence matrix
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : degree must be at least
## 1 for vertex influence matrix
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : degree must be at least
## 1 for vertex influence matrix
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : degree must be at least
## 1 for vertex influence matrixy = pred$fit
ci <- pred$se.fit * qt(0.95 / 2 + .5, pred$df)
ymin = y - ci
ymax = y + ci
loess.DF <- data.frame(x = xseq, y, ymin, ymax, se = pred$se.fit)
head(loess.DF)## x y ymin ymax se
## 1 2.000000 0.7245604 0.7183194 0.7308013 0.002852544
## 2 2.126582 0.7245711 0.7183328 0.7308095 0.002851368
## 3 2.253165 0.7246016 0.7183703 0.7308328 0.002848108
## 4 2.379747 0.7246488 0.7184283 0.7308694 0.002843217
## 5 2.506329 0.7247100 0.7185027 0.7309174 0.002837193
## 6 2.632911 0.7247824 0.7185896 0.7309751 0.002830529summary(loess.DF)## x y ymin ymax
## Min. : 2.0 Min. :0.7246 Min. :0.7183 Min. :0.7308
## 1st Qu.: 4.5 1st Qu.:0.7264 1st Qu.:0.7202 1st Qu.:0.7325
## Median : 7.0 Median :0.7279 Median :0.7221 Median :0.7337
## Mean : 7.0 Mean :0.7274 Mean :0.7216 Mean :0.7332
## 3rd Qu.: 9.5 3rd Qu.:0.7284 3rd Qu.:0.7230 3rd Qu.:0.7341
## Max. :12.0 Max. :0.7293 Max. :0.7237 Max. :0.7349
## se
## Min. :0.002205
## 1st Qu.:0.002519
## Median :0.002716
## Mean :0.002655
## 3rd Qu.:0.002802
## Max. :0.003012summary(mloess)## Call:
## loess(formula = wing2body ~ month_of_year, data = dfF, degree = 0)
##
## Number of Observations: 14
## Equivalent Number of Parameters: 2.17
## Residual Standard Error: 0.007079
## Trace of smoother matrix: 2.4 (exact)
##
## Control settings:
## span : 0.75
## degree : 0
## family : gaussian
## surface : interpolate cell = 0.2
## normalize: TRUE
## parametric: FALSE
## drop.square: FALSE# compare models
mlinear = glm(wing2body ~ month_of_year, data=df, family=gaussian)
mlinearF = glm(wing2body ~ month_of_year, data=dfF, family=gaussian)
mlinearM = glm(wing2body ~ month_of_year, data=dfM, family=gaussian)
summary(mlinear)$coeff[[8]] # significant for all (but only with K.elegans)## [1] 0.0222125summary(mlinearF)$coeff # not significant for females (or when split by host plant)## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.720672319 0.0043862061 164.304254 1.735955e-21
## month_of_year 0.000934252 0.0005570487 1.677146 1.193477e-01summary(mlinearM)$coeff # significant for males (for K.elegans, not for C.corindum)## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.726628910 0.0026151032 277.858597 3.177879e-24
## month_of_year 0.000708547 0.0003321184 2.133417 5.422375e-02scale_colour_discrete(name =“Host Plant”, labels=c(“C. corindum”, “K. elegans”)) +
Wing Morph Freq
#############wing morph plots
##Results to visualize: sex, host, months_since_start.
w_morph_summary<-aggregate(wing_morph_b~sex*pophost*month_of_year, data=raw_data, FUN=mean)
w_morph_summary$se<-aggregate(wing_morph_b~sex*pophost*month_of_year, data=raw_data, FUN=SE)$wing_morph_b
#jitter slightly
jitter = runif(n=nrow(w_morph_summary), min=-0.5, max=0.5)
w_morph_summary$dates <- w_morph_summary$month_of_year + jitter
# create groups
dd = w_morph_summary
dd$pophost[dd$pophost=="C.corindum"]<-"C. corindum"
dd$pophost[dd$pophost=="K.elegans"]<-"K. elegans"
dd$pophost = factor(dd$pophost, levels = c("K. elegans", "C. corindum") )
dd$`Host Plant` = dd$pophost
dd$sex[dd$sex=="F"]<-"Females"
dd$sex[dd$sex=="M"]<-"Males"
dd$sex = factor(dd$sex, levels = c("Males", "Females") )
dd$`Sex` = dd$sex
dfF = dd[dd$sex=="Females",]
dfM = dd[dd$sex=="Males",]p7 = ggplot() +
ggtitle("Females") + xlab("Month") + ylab("Long-Wing Morph Frequency") +
geom_vline(xintercept = xlab_dates, color="gainsboro") +
geom_smooth(data=dfF, method="lm", se=FALSE, linetype = "dashed", lwd=0.5,
mapping = aes(x = month_of_year, y = wing_morph_b, colour = `Host Plant`)) +
geom_smooth(data=dfF, method="loess",
mapping = aes(x = month_of_year, y = wing_morph_b, colour=`Host Plant`, fill=`Host Plant`)) +
geom_point(data=dfF, mapping = aes(x = month_of_year, y = wing_morph_b, colour=`Host Plant`)) +
customPlot +
scale_color_manual(values=c("C. corindum" = "turquoise3", "K. elegans" = "springgreen4")) +
scale_fill_manual(values = c("C. corindum" = "turquoise3", "K. elegans" = "green")) +
scale_x_continuous(breaks=xlab_months, labels= month_labs) +
ylim(-0.75,2.25)
p8 = ggplot() +
ggtitle("Males") + xlab("Month") + ylab("Long-Wing Morph Frequency") +
geom_vline(xintercept = xlab_dates, color="gainsboro") +
geom_smooth(data=dfM, method="lm", se=FALSE, linetype = "dashed", lwd=0.5,
mapping = aes(x = month_of_year, y = wing_morph_b, colour = `Host Plant`)) +
geom_smooth(data=dfM, method="loess",
mapping = aes(x = month_of_year, y = wing_morph_b, colour=`Host Plant`, fill=`Host Plant`)) +
geom_point(data=dfM, mapping = aes(x = month_of_year, y = wing_morph_b, colour=`Host Plant`)) +
customPlot +
scale_color_manual(values=c("C. corindum" = "turquoise3", "K. elegans" = "springgreen4")) +
scale_fill_manual(values = c("C. corindum" = "turquoise3", "K. elegans" = "green")) +
scale_x_continuous(breaks=xlab_months, labels= month_labs) +
ylim(-0.75,2.25)
grid.arrange(p7,p8, ncol=2)## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
extrafigure = grid.arrange(p3,p4,p7,p8, ncol=2)## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
#ggsave("extra_figure.pdf", plot=extrafigure)customPlot = list( theme_classic(),
theme(axis.text=element_text(size=13),
axis.title=element_text(size=16),
plot.title=element_text(size=20),),
theme(legend.position = c(0.23, 0.95)),
theme(legend.title = element_text(size=14, face="italic"),
legend.text = element_text(size = 13, face="italic"))
)#########
p9 = ggplot() +
ggtitle("A") + xlab("Month") + ylab("Long-Wing Morph Frequency") +
geom_vline(xintercept = xlab_dates, color="gainsboro") +
geom_smooth(data=dd, method="lm", se=FALSE, linetype = "dashed", lwd=0.5,
mapping = aes(x = month_of_year, y = wing_morph_b), colour="black") +
geom_smooth(data=dd, method="loess",
mapping = aes(x = month_of_year, y = wing_morph_b, colour=`Host Plant`, fill=`Host Plant`)) +
geom_point(data=dd, mapping = aes(x = month_of_year, y = wing_morph_b, colour=`Host Plant`)) +
customPlot +
scale_color_manual(values=c("C. corindum" = "turquoise3", "K. elegans" = "springgreen4")) +
scale_fill_manual(values = c("C. corindum" = "turquoise3", "K. elegans" = "green")) +
scale_x_continuous(breaks=xlab_months, labels= month_labs) + ylim(0.3,1.2)
p10 = ggplot() +
ggtitle("B") + xlab("Month") + ylab("Long-Wing Morph Frequency") +
geom_vline(xintercept = xlab_dates, color="gainsboro") +
geom_smooth(data=dd, method="lm", se=FALSE, linetype = "dashed", lwd=0.5,
mapping = aes(x = month_of_year, y = wing_morph_b), colour="black") +
geom_smooth(data=dd, method="loess",
mapping = aes(x = month_of_year, y = wing_morph_b, colour=Sex, fill=Sex)) +
geom_point(data=dd, mapping = aes(x = month_of_year, y = wing_morph_b, colour=Sex)) +
customPlot +
scale_color_manual(values=c("Females" = "brown3", "Males" = "black")) +
scale_fill_manual(values = c("Females" = "brown1", "Males" = "sienna4")) +
scale_x_continuous(breaks=xlab_months, labels= month_labs) + theme(axis.title.y = element_blank()) + ylim(0.3,1.2)
# extract the LOESS
alpha = paste("alpha[loess]==", ggplot_build(p9)$data[[2]]$alpha[1])## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'degree="lambda[loess]==0"
mlinear = glm(wing_morph_b ~ month_of_year, data=dd, family=gaussian)
round(summary(mlinear)$coeff[[8]],3) # significant for all (but only with K.elegans)## [1] 0.487pvalue = paste0("italic(p)[glm]==",round(summary(mlinear)$coeff[[8]],2))
p9 = p9 + annotate(geom="text", x=7.4, y=1.0, label=alpha, color="black", parse=TRUE) +
annotate(geom="text", x=7.4, y=1.04, label=degree, color="black", parse=TRUE) +
annotate(geom="text", x=5, y=0.718, label=pvalue, color="black", parse=TRUE)
p10 = p10 + annotate(geom="text", x=7.4, y=1.1, label=alpha, color="black", parse=TRUE) +
annotate(geom="text", x=7.4, y=1.15, label=degree, color="black",parse=TRUE) +
annotate(geom="text", x=5, y=1, label=pvalue, color="black", parse=TRUE)
grid.arrange(p9,p10, ncol=2)## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
threshold = sort(unique(raw_data$dates))[5]
df$date_b = "2013-2015"
df$date_b[df$dates >= threshold] = "2016-2020"
df$date_b = as.factor(df$date_b)wm_table<-aggregate(wing_morph_b~dates, data=raw_data, FUN=mean)
wm_table$date_b = "2013-2015"
wm_table$date_b[wm_table$dates >= threshold] = "2016-2020"
wm_table$date_b = as.factor(wm_table$date_b)
last_4_years = wm_table[wm_table$date_b=="2016-2020",]
first_3_years = wm_table[wm_table$date_b=="2013-2015",]
xlab_years = sort(unique(raw_data$dates))p11 = ggplot() +
ggtitle("C") + xlab("Year") + ylab("Long-Wing Morph Frequency") +
geom_vline(xintercept = xlab_years, color="gainsboro") +
geom_smooth(data=wm_table, method="lm", se=FALSE, linetype = "dashed",
mapping = aes(x = dates, y = wing_morph_b), colour="black", lwd=0.5) +
geom_smooth(data=wm_table, method="loess",
mapping = aes(x = dates, y = wing_morph_b), colour="black") +
geom_point(data=wm_table, mapping = aes(x = dates, y = wing_morph_b))
#ylim(0.71, 0.75) +
#+ stat_smooth(data=first_3_years, mapping = aes(x = dates, y = wing_morph_b), method = "lm", formula = y ~ x, se = FALSE,
# colour="blue") # linetype = "longdash",
p11 = p11 + guides(fill = guide_legend(reverse = TRUE)) + theme_classic() +
theme(axis.text=element_text(size=13),
axis.title=element_text(size=16),
plot.title=element_text(size=20)) + guides(color = FALSE) + scale_linetype(guide = FALSE) +
theme(legend.key = element_rect(fill = "white", color = NA),
legend.key.size = unit(0.5, "cm"),
legend.key.width = unit(0.5,"cm")) +
scale_fill_manual(values = "blue", labels=c("C. corindum", "K. elegans")) +
labs(fill = "Host Plant") + theme(legend.title = element_text(size = 14),
legend.text = element_text(size = 13, face="italic")) +
scale_color_manual(values="blue") + theme(legend.position = c(0.2, 0.88)) + theme(legend.title = element_blank()) + ylim(0.3,1.2)
alpha = paste("alpha[loess]==", ggplot_build(p0)$data[[2]]$alpha[1])## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'degree="lambda[loess]==0"
mlinear = glm(wing_morph_b ~ dates, data=dd, family=gaussian)
pvalue = paste0("italic(p)[glm]==",round(summary(mlinear)$coeff[[8]],2))
p11 = p11 + annotate(geom="text", x=unique(raw_data$dates)[7], y=1, label=alpha, color="black", parse=TRUE, size=6) +
annotate(geom="text", x=unique(raw_data$dates)[7], y=1.05, label=degree, color="black",parse=TRUE, size=6) +
annotate(geom="text", x=unique(raw_data$dates)[5], y=0.59, label=pvalue, color="black", parse=TRUE, size=6)
# first 3 years
#mlinear = glm(wing_morph_b ~ dates, data=first_3_years, family=gaussian)
#pvalue = paste0("italic(p)[glm]==",round(summary(mlinear)$coeff[[8]],2))
#p11 = p11 + annotate(geom="text", x=unique(raw_data$dates)[3], y=1, label=pvalue, color="blue", parse=TRUE, size=6)
p11## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
figure2 = ggdraw() +
draw_plot(p9, 0, .5, .5, .5) +
draw_plot(p10, 0.5, .5, .5, .5) +
draw_plot(p11, 0, 0, 1, .5)## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'figure2
#ggsave("wmorph_over_time.pdf", plot=figure2)MLC: Something is going wrong here calculating variance for a binomial response - a frequency or probability plot should have a y-axis with a strict minimum of 0 and maximum of 1, but that is not happening with the variance plotted here. So, the estimated variance is probably not accurate, and the expanded y-axis range makes it difficult to see changes that might be visible along the month axis. I have historically used binom.confint() in the binom() library to calculate binomial confidence intervals; I don’t know what the smoothing options are, but I do think we need a different function for this one. In principal, I like this idea.
MLC: Actually, this is cuing me into something being off about the way loess is visualizing variance in general - the error bars on many of these plots are way wider (x10?) than the standard error on the data itself. I think it may be related to using a summary file instead of a complete data file.
AB:
The prediction interval predicts in what range a future individual observation will fall, while a confidence interval shows the likely range of values associated with some statistical parameter of the data, such as the population mean.
What is being plotted is a CI but a CI will converge to a single value as the sample size increases. That is why the shaded regions change so much as you keep adding more data points.
What we want then is a prediction interval. Even though the “y” column is the predicted value, the “ymin” and “ymax” are the lower and upper pointwise confidence interval around the mean, note the predicted interval.
library(msir)## Package 'msir' version 1.3.3## Type 'citation("msir")' for citing this R package in publications.data(cars)
# Calculates and plots a 1.96 * SD prediction band, that is,
# a 95% prediction band
l <- loess.sd(cars, nsigma = 1.96)
plot(cars, main = "loess.sd(cars)", col="red", pch=19)
lines(l$x, l$y)
lines(l$x, l$upper, lty=2)
lines(l$x, l$lower, lty=2)