--- title: "EcoStage_mixedmodeling_Apr2020" author: "Alycia Lackey" date: "April 20, 2020" output: html_document --- ======================================== OFFSPRING GROWTH & SIZE, PRE & POST-DIET ======================================== ```{r} # PREDIET # import data library(readxl) prediet_data <- read_excel("Growth_data_ARL.xlsx", na = ".", col_types = c("text", "text", "text", "text", "text", "text", "text", "numeric", "numeric", "numeric")) View(prediet_data) # visualize data, histogram library(ggplot2) ggplot(prediet_data, aes(x=prediet_growth)) + geom_histogram() ggplot(prediet_data, aes(x=prediet_growth)) + geom_histogram() + facet_wrap(~species*dadcolor) # check normality, NS p-value qqnorm(prediet_data$prediet_growth) qqline(prediet_data$prediet_growth) shapiro.test(prediet_data$prediet_growth) # model testing library(lme4) prediet_growth_lmer <- lmer(prediet_growth ~ species * dadcolor + (1|dadID) + (1|momID), data = prediet_data, family = "gaussian") summary(prediet_growth_lmer) library(lmerTest) prediet_g <- lmer(prediet_growth ~ species * dadcolor + (1|dadID) + (1|momID), data = prediet_data) summary(prediet_g) ``` ```{r} # postDIET # import data library(readxl) postdiet_data <- read_excel("Growth_data_ARL.xlsx", na = ".", sheet = "Post_diet", col_types = c("text", "text", "text", "text", "text", "text", "text", "text", "numeric", "numeric", "numeric")) View(postdiet_data) # histogram library(ggplot2) ggplot(postdiet_data, aes(x=postdiet_growth)) + geom_histogram() ggplot(postdiet_data, aes(x=postdiet_growth)) + geom_histogram() + facet_wrap(~species*dadcolor*diet) # check normality, NS p-value qqnorm(postdiet_data$postdiet_growth) qqline(postdiet_data$postdiet_growth) shapiro.test(postdiet_data$postdiet_growth) # model testing library(lme4) postdiet_growth_glmer <- glmer(postdiet_growth ~ species * dadcolor * diet + (1|dadID) + (1|momID), data = postdiet_data, family = "gaussian") summary(postdiet_growth_glmer) library(lmerTest) postdiet_g <- lmer(postdiet_growth ~ species * dadcolor * diet + (1|dadID) + (1|momID), data = postdiet_data) summary(postdiet_g) ``` ```{r} # Size at Maturity (48 week size) # histogram library(ggplot2) ggplot(postdiet_data, aes(x=avglen_48wks)) + geom_histogram() ggplot(postdiet_data, aes(x=avglen_48wks)) + geom_histogram() + facet_wrap(~species*dadcolor*diet) # check normality, NS p-value qqnorm(postdiet_data$avglen_48wks) qqline(postdiet_data$avglen_48wks) shapiro.test(postdiet_data$avglen_48wks) # model testing library(lme4) avgleng48wks_glmer <- glmer(avglen_48wks ~ species * dadcolor * diet + (1|dadID) + (1|momID), data = postdiet_data, family = "gaussian") summary(avgleng48wks_glmer) library(lmerTest) avgleng48wks_lmertest <- lmer(avglen_48wks ~ species * dadcolor * diet + (1|dadID) + (1|momID), data = postdiet_data) summary(avgleng48wks_lmertest) ``` ======================================== SON FITNESS (BODY SIZE, ATTRACTIVENESS, FEMALE RESPONSE/PREFERENCE) ======================================== ```{r} # Import data library(readxl) msuccess <- read_excel("msuccess06new_RT_ARL_10.30.19.xlsx", sheet = "msuccessRtrim", col_types = c("text", "text", "text", "numeric", "text", "text", "text", "numeric", "numeric", "numeric", "numeric", "numeric", "numeric", "numeric", "numeric", "numeric", "numeric", "numeric", "numeric", "numeric", "numeric", "numeric", "numeric", "numeric"), na = ".") View(msuccess) attach(msuccess) ``` Visualize son metrics of mating success ```{r} hist(msuccess$pref) hist(msuccess$folperleadbin) hist(msuccess$courtshipvigor) hist(msuccess$logcourtshipv) #normal distribution of log-transformed & qqnorm test cor.test(avgredindex, logcourtshipv, use = "complete.obs") library("ggpubr") ggscatter(msuccess, x = "avgredindex", y = "logcourtshipv", add = "reg.line", conf.int = TRUE, cor.coef = TRUE, cor.method = "pearson", xlab = "red", ylab = "courtvigor") ``` Normality testing ```{r} # log transf of courtship vigor qqnorm(msuccess$logcourtshipv) qqline(msuccess$logcourtshipv) shapiro.test(msuccess$logcourtshipv) # pref score, histogram is bell shaped, but qqplot is non-normal, use poisson qqnorm(msuccess$pref) qqline(msuccess$pref) shapiro.test(msuccess$pref) ``` ---Male mating success variables--- folperleadbin (binomial) logcourtshipvigor (normal) preference score (poisson) random effect = trial nested within male ```{r} library(lme4) library(lmerTest) library(car) library(ggplot2) #logcourtshipv logcourt_lmm <- lmer(logcourtshipv ~ sp * dadc * diet + (trial|male), msuccess) summary(logcourt_lmm) ggplot(msuccess, aes(x=sp, y=logcourtshipv)) + geom_boxplot() ggplot(msuccess, aes(x=sp, y=logcourtshipv, fill=dadc)) + geom_boxplot() library(lmerTest) court_fit <- lmer(logcourtshipv ~ sp * dadc * diet + (trial|male), msuccess) anova(court_fit, type = "III", dff = "Satterthwaite") logcourt_lsmeans <- lsmeans(court_fit, pairwise ~ sp:dadc, adjust = "none") #adjusting p-value pcourt <- c(0.49, 0.0600, 0.4284, 0.0005) p.adjust(pcourt, method = "fdr", n = 4) ``` ```{r} # follow per lead (binary) follead_glmm <- glmer(folperleadbin ~ sp * dadc * diet + (trial|male), data = msuccess, family = "binomial") summary(follead_glmm) folleadmtraits_glmm <- glmer(folperleadbin ~ avgredindex + logcourtshipv + (trial|male), data = msuccess, family = "binomial") summary(folleadmtraits_glmm) folleadmtraits_glm <- glm(folperleadbin ~ avgredindex + logcourtshipv, data = msuccess, family = "binomial") summary(folleadmtraits_glm) ``` ```{r} library(lme4) # preference score pref_rep <- glmer(pref ~ sp * dadc * diet + (trial|male), data = msuccess, family = "poisson") summary(pref_rep) # pref score predicted by male traits pref_repmtraits <- glmer(pref ~ avgredindex + logcourtshipv + (trial|male), data = msuccess, family = "poisson") summary(pref_repmtraits) ``` ---Male attractiveness--- ```{r} library(readxl) attract <- read_excel("C:/Users/alyci/OneDrive/Desktop/EcoStage/attract_11_08_19.xlsx", col_types = c("text", "text", "text", "text", "text", "text", "text", "numeric", "numeric", "numeric", "text")) View(attract) hist(attract$redindex) hist(attract$avglen) hist(attract$mwt) qqnorm(attract$redindex) qqline(attract$redindex) shapiro.test(attract$redindex) qqnorm(attract$avglen) qqline(attract$avglen) shapiro.test(attract$avglen) qqnorm(attract$mwt) qqline(attract$mwt) shapiro.test(attract$mwt) ``` ---------------------------------------------------------------------------- REDINDEX FROM ATTRACT DATA SET (all males, incl those without mating trials) ---------------------------------------------------------------------------- ```{r} library(lme4) library(car) library(ggplot2) library(emmeans) # redinex, random effect of cross ID - controls for mom & dad red_attract_lmm <- lmer(attract$redindex ~ attract$species * attract$dad_color * attract$diet + (1|attract$xid)) summary(red_attract_lmm) # boxplot ggplot(attract, aes(x=diet, y=redindex, fill=dad_color)) + geom_boxplot() + facet_wrap(~species) # CONTRASTS # make categorical variables factors attract$species.factor <- factor(attract$species) is.factor(attract$species.factor) attract$dadcol.factor <- factor(attract$dad_color) is.factor(attract$dadcol.factor) attract$diet.factor <- factor(attract$diet) is.factor(attract$diet.factor) levels(attract$species.factor) levels(attract$dadcol.factor) levels(attract$diet.factor) red_attract_model <- lmer(attract$redindex ~ attract$species.factor * attract$dadcol.factor * attract$diet.factor + (1|attract$xid)) summary(red_attract_model) attach(attract) red_attract_model1 <- lmer(redindex ~ species.factor * dadcol.factor * diet.factor + (1|xid)) summary(red_attract_model1) red_lsmeans <- lsmeans(red_attract_model1, pairwise ~ species.factor:dadcol.factor:diet.factor, adjust = "none") sp.red <- c("PB", "PL", "PB", "PL", "PB", "PL", "PB", "PL") dadc.red <- c("D", "D", "H", "H", "D", "D", "H", "H") diet.red <- c("i", "i", "i", "i", "p", "p", "p", "p") lsmean.red <- c(1.66, 3.43, 2.51, 2.33, 1.44, 1.85, 1.38, 1.83) SE.red <- c(0.273, 0.434, 0.282, 0.341, 0.310, 0.438, 0.254, 0.316) df.red <- c(33, 69, 51, 70, 50, 103, 34, 61) LCL.red <- c(1.102, 2.564, 1.943, 1.646, 0.813, 0.981, 0.866, 1.196) UCL.red <- c(2.21, 4.29, 3.08, 3.01, 2.06, 2.72, 1.90, 2.46) df.lsmeans.red <- data.frame(sp.red, dadc.red, diet.red, lsmean.red, SE.red, df.red, LCL.red, UCL.red) View(df.lsmeans.red) sp.labs <- c(PB = "Benthic", PL = "Limnetic") ggplot(data=df.lsmeans.red, aes(x=diet.red, y=lsmean.red, fill=dadc.red)) + geom_point(color = c("#C2160A", "#F0B8B4", "#C2160A", "#F0B8B4", "#C2160A", "#F0B8B4", "#C2160A", "#F0B8B4"), position = position_dodge(width = 0.4), size = 5) + geom_errorbar(aes(ymin=LCL.red, ymax=UCL.red), color = c("#C2160A", "#F0B8B4", "#C2160A", "#F0B8B4", "#C2160A", "#F0B8B4", "#C2160A", "#F0B8B4"), width = 0.4, position = "dodge") + facet_grid(~sp.red, labeller = labeller(sp.red = sp.labs)) + xlab("diet") + ylab ("mean nuptial color") + theme(panel.grid.major = element_blank(), panel.background = element_rect(fill = "white"), axis.line = element_line(size = 0.5), axis.title = element_text(size = 18, face = "bold"), axis.title.y = element_text(margin = margin(t = 0, r = 10, b = 0, l = 0, unit = "pt")), axis.title.x = element_text(margin = margin(t = 0, r = 0, b = 10, l = 0, unit = "pt")), axis.text = element_text(size = 16, colour = "black"), legend.text = element_text(size = 13), legend.title = element_blank(), legend.spacing.x = unit(0.1, "cm"), strip.text.x = element_text(size = 15, face = "bold")) #adjusting p-value p <- c(0.0011, 0.6820, 0.4439, 0.2768, 0.0358, 0.08919, 0.0498, 0.9669, 0.4804, 0.0001, 0.0037, 0.2119) p.adjust(p, method = "fdr", n = 12) ``` ```{r} library(lme4) library(car) library(ggplot2) #length w random effect of cross ID length_attract_lmm <- lmer(attract$avglen ~ attract$species * attract$dad_color * attract$diet + (1|attract$xid)) summary(length_attract_lmm) length_lsmeans <- lsmeans(length_attract_lmm, pairwise ~ species, adjust = "none") ggplot(attract, aes(x=species, y=avglen)) + geom_boxplot() ``` ```{r} library(lme4) library(car) library(ggplot2) #weight, random effect of cross ID wt_attract_lmm <- lmer(attract$mwt ~ attract$species * attract$dad_color * attract$diet + (1|attract$xid)) summary(wt_attract_lmm) wt_lsmeans <- lsmeans(wt_attract_lmm, pairwise ~ species, adjust = "none") ggplot(attract, aes(x=species, y=mwt)) + geom_boxplot() ``` ```{r} cor.test(attract$avglen, attract$mwt) library("ggpubr") ggscatter(attract, x = "avglen", y = "mwt", add = "reg.line", conf.int = TRUE, cor.coef = TRUE, cor.method = "pearson", xlab = "length", ylab = "wt") library(lmf) extract <- lm.extract(formula = attract$mwt ~ attract$redindex, data = attract) str(extract) lenwtresid <- extract$res lenwtresid attract$lenwtresid <- lenwtresid str(attract) hist(attract$lenwtresid) #length-weight residuals, random effect of cross ID lwresid_attract_lmm <- lmer(attract$lenwtresid ~ attract$species * attract$dad_color * attract$diet + (1|attract$xid)) summary(lwresid_attract_lmm) lwresid_lsmeans <- lsmeans(lwresid_attract_lmm, pairwise ~ species, adjust = "none") ggplot(attract, aes(x=species, y=lenwtresid)) + geom_boxplot() ggplot(attract, aes(x=dad_color, y=lenwtresid)) + geom_boxplot() ```