##loading packages## library(EnvStats) library(ggplot2) library(ggpubr) library(ggfortify) library(dplyr) library(vegan) library(olsrr) library(MuMIn) library(MASS) ##for glm## library(pscl) #Importing data and making row 1 the row names data <- read.csv("cleandata.csv",row.names = 1) ##data proccessing for sperm and polyp counts, plus normalization for symbiont data## #averaging sperm into new column data$spermavg <- rowMeans(subset(data, select = c(spermcount1, spermcount2, spermcount3))) data$spermtot = data$spermavg*1000000 #averaging polyp counts data$polypavg <- rowMeans(subset(data, select = c(polypcount1, polypcount2, polypcount3))) #deleting count columns data <- data[-c(10:15,18:21)] #standardizing sperm/size data$spermsizestand <- data$spermtot/data$size #standardizing sperm/polyp data$spermpolypstand<- data$spermtot/data$polypavg ##outlier tests and removals## rosnerTest(data$cat_stand) ##none## rosnerTest(data$pox_stand) ##one outlier- RI55## rosnerTest(data$ppo_stand) ##none## rosnerTest(data$dt) ##none## rosnerTest(data$mel_stand) ##one outlier- RI51## rosnerTest(data$lipid_stand) ##one outlier- RI53## rosnerTest(data$carb_stand) ##no outlier## rosnerTest(data$symbiont_stand) ##no outlier## rosnerTest(data$spermpolypstand) ##no outlier## #Getting rid of outliers in full set## data_noout <- data[!(row.names(data) %in% c("RI53", "RI51", "RI55")),] ##tests for normality and transformations## shapiro.test(data_noout$cat_stand) shapiro.test(data_noout$pox_stand) shapiro.test(data_noout$ppo_stand) ##not normal## shapiro.test(data_noout$dt) shapiro.test(data_noout$mel_stand) shapiro.test(data_noout$carb_stand) shapiro.test(data_noout$lipid_stand) shapiro.test(data_noout$symbiont_stand) ##not normal.. and will never be, but not an issue as it's always a predictor## shapiro.test(data_noout$spermpolypstand) ##not normal## ##plot non-normal data## hist(data_noout$ppo_stand, main = "PPO", breaks = 10) ##transform data_noout$ppo_stand.norm <- log(data_noout$ppo_stand) ##check transformed## shapiro.test(data_noout$ppo_stand.norm) ##fixed## hist(data_noout$spermpolypstand, main = "Sperm", breaks = 10) ##zero inflated## data_noout$spermpolypstand.norm <- log(data_noout$spermpolypstand) ##check transformed## shapiro.test(data_noout$spermpolypstand.norm) ##fixed## ##hypothesis testing time## ##start by confirming that white and brown corals have different symbiont densities## sym_test=data[,c(1,10)] wilcox.test(data$symbiont_stand~data$type) ####SIGNIFICANT, p < 0.001 ##make a figure## sym_test %>% ggplot( aes(x=factor(type, level = c('White', 'Brown')), y=symbiont_stand, fill=type)) + geom_boxplot() + geom_jitter(color="black", size=2, alpha=0.9) + theme_classic()+ scale_fill_manual(values = c('Brown' = "burlywood4", 'White' ="lightgrey")) + theme( legend.position="none") + xlab("Symbiotic State") + ylab("Symbiont Density") ##check for effects of symbiont density on lipids/carbs## sym_lip_cont=data[,c(8,10)] sym_carb_cont=data[,c(9,10)] sym_lip_bi=data[,c(1,8)] sym_carb_bi=data[,c(1,9)] sym_lip_cont=na.omit(sym_lip_cont) sym_lip_bi=na.omit(sym_lip_bi) ##test for outliers/normality## rosnerTest(sym_lip_cont$lipid_stand) rosnerTest(sym_lip_bi$lipid_stand) sym_lip_cont <- sym_lip_cont[!(row.names(sym_lip_cont) %in% c("RI53")),] sym_lip_bi <- sym_lip_bi[!(row.names(sym_lip_bi) %in% c("RI53")),] rosnerTest(sym_carb_cont$carb_stand) rosnerTest(sym_carb_bi$carb_stand) shapiro.test(sym_lip_cont$lipid_stand) shapiro.test(sym_carb_cont$carb_stand) ##and check whether lipids/carbs were significantly affected by symbiont density## summary(lm(lipid_stand ~ symbiont_stand, data=sym_lip_cont)) summary(lm(carb_stand ~ symbiont_stand, data=sym_carb_cont)) summary(lm(lipid_stand ~ type, data=sym_lip_bi)) summary(lm(carb_stand ~ type, data=sym_carb_bi)) ####not significant either way## ##finally test for symbiont density and energetics as predictors of sperm output## sperm_out=data[,c(1,8:10,15)] sperm_out=na.omit(sperm_out) ##test for outliers/normality rosnerTest(sperm_out$spermpolypstand)#none shapiro.test(sperm_out$spermpolypstand) ##non-normal# sperm_out$spermpolypstand.log=log(sperm_out$spermpolypstand) ##fixes it## ##modeling## ##first model with symbiont density continuous## sperm_out_continuous=sperm_out[,c(2:4,6)] ##set up for dredging sperm_out_global_model<- glm(spermpolypstand.log~ . + .^3, data =sperm_out_continuous, na.action = "na.fail") #dredge all models sperm_out.all.dredge <- dredge(sperm_out_global_model, beta = c("none"), evaluate = TRUE, rank = "AICc", cluster = "NULL") #view sperm_outalase model output sperm_out.all.dredge summary(sperm_out.all.dredge) #best model is carb## summary(glm(spermpolypstand.log ~ carb_stand, data=sperm_out_continuous)) # Estimate Std. Error t value Pr(>|t|) #(Intercept) 12.2660 0.6867 17.861 1.88e-12 *** #carb_stand 224.1498 84.4700 2.654 0.0167 * #--- #Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 ##carb sig## ##then model with categorical symbiont density## sperm_out_bi=sperm_out[,c(1:3,6)] sperm_out_bi = na.omit(sperm_out_bi) ##set up for dredging sperm_out_global_model<- glm(spermpolypstand.log ~ . + .^3, data =sperm_out_bi, na.action = "na.fail") #dredge all models sperm_out.all.dredge <- dredge(sperm_out_global_model, beta = c("none"), evaluate = TRUE, rank = "AICc", cluster = "NULL") #view sperm_outalase model output sperm_out.all.dredge summary(sperm_out.all.dredge) #average models between 0-2 delta AICc sperm_out.all.model.avg <- model.avg(sperm_out.all.dredge, subset = delta <= 2) sperm_out.all.model.avg summary(sperm_out.all.model.avg) #Model-averaged coefficients: #(full average) #Estimate Std. Error Adjusted SE z value Pr(>|z|) #(Intercept) 12.1470 0.7353 0.7893 15.389 <2e-16 *** #carb_stand 218.8971 85.0768 91.6521 2.388 0.0169 * #typeWhite 0.2293 0.5152 0.5375 0.427 0.6697 #(conditional average) #Estimate Std. Error Adjusted SE z value Pr(>|z|) #(Intercept) 12.1470 0.7353 0.7893 15.389 <2e-16 *** #carb_stand 218.8971 85.0768 91.6521 2.388 0.0169 * #typeWhite 0.7847 0.6875 0.7436 1.055 0.2913 #--- #Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 ##still just carb sig## ##plot it## ggplot(sperm_out, aes(spermpolypstand.log, carb_stand)) + geom_smooth(method = "lm", color = "black", lty = 1) + geom_point(aes(fill = type), pch = 21, color = "black", size = 3) + scale_fill_manual(values = c('Brown' = "burlywood4", 'White' = "white")) + theme_bw() + theme(panel.grid = element_blank()) + ylab("Carbohdyrate Concentration\n(mg lipid/mg tissue)") + xlab("Sperm Per Polyp\n(log normalized)") ##let's try a PCA## ##remove rows with missing data## data_noout = na.omit(data_noout) #filter data frames to just numeric predictor and response variables and scale and center the variables names(data_noout) numeric_data <- data_noout[,c(3:10,15)] numeric_data_scaled <- as.data.frame(scale(numeric_data, center = TRUE, scale = TRUE)) #check for proper centering (values should be ~0) apply(numeric_data_scaled, 2, FUN = mean) #filter scaled data into predictor and response variables names(numeric_data_scaled) predictorVariables <- numeric_data_scaled[,c(6:9)] responseVariables <- numeric_data_scaled[,c(1:5)] #run PCA PCA <- prcomp(responseVariables, scale = TRUE, center = TRUE) summary(PCA) #plot PCA based on immune parameters autoplot(PCA, data = data_noout, colour = "type", loadings = TRUE, loadings.colour = 'black', loadings.label = TRUE, loadings.label.colour= "black", loadings.label.size = 3) + theme_bw() + theme(panel.grid = element_blank()) + stat_ellipse(aes (color = type))+ scale_color_manual(values = c('Brown' = "burlywood4", 'White' ="lightgrey")) ##now do an RDA #add categorical variables to predictor dataframe for RDA analysis (make sure rows are in same order to use cbind) #predictorVariables <- cbind(predictorVariables, data_noout$type) #fix column header names names(predictorVariables) #run RDA global model RDA_global <- rda(responseVariables ~ ., data = predictorVariables) summary(RDA_global) plot(RDA_global, display = c("sp", "lc", "cn")) #test significance of global model anova(RDA_global, permutations = 999) #global model significant, p = 0.027 #test significance by axis anova(RDA_global, by = "axis", permutations = 999) #no RDA axes are significant #test significance by margin anova(RDA_global, by = "margin", permutations = 999) # carb is significant (0.042) p-values in parentheses #test significance by term anova(RDA_global, by = "term", permutations = 999) #spermpolypstand is significant (0.039), p-values in parentheses #check for colinearity vif.cca(RDA_global) #use forward calling to identify significant predictor variables RDA_forwardSel <- ordiR2step(rda(responseVariables ~ 1, data = predictorVariables), scope = formula(RDA_global), direction = "forward", R2scope = TRUE, permutations = 999) ##nothing is significant.. do not continue with RDA## ##now we'll try the correlation matrix## ##for this we'll leave in the outliers to start## corrdata=data_noout corrdata=corrdata[,c(3:4,6,16,7:10,17)] library(corrplot) ##function for p values## cor.mtest <- function(mat, ...) { mat <- as.matrix(mat) n <- ncol(mat) p.mat<- matrix(NA, n, n) diag(p.mat) <- 0 for (i in 1:(n - 1)) { for (j in (i + 1):n) { tmp <- cor.test(mat[, i], mat[, j], ...) p.mat[i, j] <- p.mat[j, i] <- tmp$p.value } } colnames(p.mat) <- rownames(p.mat) <- colnames(mat) p.mat } ##run correlations and p values M<-cor(corrdata) p.mat <- cor.mtest(corrdata) ##plot## col <- colorRampPalette(c("#4477AA","#77AADD","#FFFFFF","#EE9988", "#BB4444")) corrplot(M, method="color", col=col(200), type="lower", addgrid.col = TRUE, addCoef.col = "black", # Add coefficient of correlation tl.col="black", tl.srt=45, #Text label color and rotation # Combine with significance p.mat = p.mat, sig.level = 0.05, insig = "blank", # hide correlation coefficient on the principal diagonal diag=TRUE ) ##sperm is + corr with ppo, carb is + corr with dt and lipid, sym is + corr with mel## ##make histograms to add in to that plot## cat=ggplot(corrdata, aes(x=cat_stand)) + geom_histogram(binwidth=50) + theme_classic() + theme(axis.title.x=element_blank(),axis.title.y=element_blank()) + scale_x_continuous(expand = c(.01, 0)) + scale_y_continuous(expand = c(0, 0)) pox=ggplot(corrdata, aes(x=pox_stand)) + geom_histogram(binwidth=.2) + theme_classic() + theme(axis.title.x=element_blank(),axis.title.y=element_blank()) + scale_x_continuous(expand = c(.01, 0)) + scale_y_continuous(expand = c(0, 0)) dt=ggplot(corrdata, aes(x=dt)) + geom_histogram(binwidth=.05) + theme_classic() + theme(axis.title.x=element_blank(),axis.title.y=element_blank()) + scale_x_continuous(expand = c(.01, 0)) + scale_y_continuous(expand = c(0, 0)) ppo=ggplot(corrdata, aes(x=ppo_stand.norm)) + geom_histogram(binwidth=.1) + theme_classic() + theme(axis.title.x=element_blank(),axis.title.y=element_blank()) + scale_x_continuous(expand = c(.01, 0)) + scale_y_continuous(expand = c(0, 0)) mel=ggplot(corrdata, aes(x=mel_stand)) + geom_histogram(binwidth=.000005) + theme_classic() + theme(axis.title.x=element_blank(),axis.title.y=element_blank()) + scale_x_continuous(expand = c(.01, 0)) + scale_y_continuous(expand = c(0, 0)) lip=ggplot(corrdata, aes(x=lipid_stand)) + geom_histogram(binwidth=.0005) + theme_classic() + theme(axis.title.x=element_blank(),axis.title.y=element_blank()) + scale_x_continuous(expand = c(.01, 0)) + scale_y_continuous(expand = c(0, 0)) carb=ggplot(corrdata, aes(x=carb_stand)) + geom_histogram(binwidth=.0006) + theme_classic() + theme(axis.title.x=element_blank(),axis.title.y=element_blank()) + scale_x_continuous(expand = c(.01, 0)) + scale_y_continuous(expand = c(0, 0)) sym=ggplot(corrdata, aes(x=symbiont_stand)) + geom_histogram(binwidth=13250) + theme_classic() + theme(axis.title.x=element_blank(),axis.title.y=element_blank()) + scale_x_continuous(expand = c(.01, 0)) + scale_y_continuous(expand = c(0, 0)) sperm=ggplot(corrdata, aes(x=spermpolypstand.norm)) + geom_histogram(binwidth=.25) + theme_classic() + theme(axis.title.x=element_blank(),axis.title.y=element_blank()) + scale_x_continuous(expand = c(.01, 0)) + scale_y_continuous(expand = c(0, 0)) ggarrange(cat, pox, dt, ppo, mel, lip, carb, sym, sperm, nrow = 9, ncol =1, align = "hv", common.legend = TRUE, legend="right") ##the rest is for looking at effects of factors on immunity (symbiont density, energetics, sperm output) ##MANOVA of differences in immune parameters between whites and browns## res.man <- manova(cbind(cat_stand,pox_stand,dt,mel_stand,ppo_stand.norm) ~ type, data = data_noout) summary(res.man) ##not significant-none of our multivariate stats are significant (RDA or MANOVA), so we'll move on to univariate stats## ##we'll do univariate stats for each immune assay with model dredging/averaging## ##start with catalase## cat=data[,c(1,3,8:10,15)] ##test for outliers and normality## rosnerTest(cat$cat_stand) ##none## shapiro.test(cat$cat_stand) ##normal## ##first model with symbiont density continuous## cat_continuous=cat[,c(2:6)] cat_continuous = na.omit(cat_continuous) ##set up for dredging cat_global_model<- glm(cat_stand ~ . + .^4, data =cat_continuous, na.action = "na.fail") #dredge all models cat.all.dredge <- dredge(cat_global_model, beta = c("none"), evaluate = TRUE, rank = "AICc", cluster = "NULL") #view catalase model output cat.all.dredge summary(cat.all.dredge) #average models between 0-2 delta AICc cat.all.model.avg <- model.avg(cat.all.dredge, subset = delta <= 2) cat.all.model.avg summary(cat.all.model.avg) #Model-averaged coefficients: #(full average) #Estimate Std. Error Adjusted SE z value Pr(>|z|) #(Intercept) 1386.49 134.53 141.53 9.797 <2e-16 *** #spermpolypstand -11.04 20.17 20.86 0.529 0.597 #carb_stand 9828.53 18195.82 18819.96 0.522 0.602 #(conditional average) #Estimate Std. Error Adjusted SE z value Pr(>|z|) #(Intercept) 1386.49 134.53 141.53 9.797 <2e-16 *** #spermpolypstand -29.73 23.24 24.83 1.198 0.231 #carb_stand 26781.20 21169.55 22607.55 1.185 0.236 #--- #Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 ##then model with categorical symbiont density## cat_bi=cat[,c(1:4,6)] cat_bi = na.omit(cat_bi) ##set up for dredging cat_global_model<- glm(cat_stand ~ . + .^4, data =cat_bi, na.action = "na.fail") #dredge all models cat.all.dredge <- dredge(cat_global_model, beta = c("none"), evaluate = TRUE, rank = "AICc", cluster = "NULL") #view catalase model output cat.all.dredge summary(cat.all.dredge) #average models between 0-2 delta AICc cat.all.model.avg <- model.avg(cat.all.dredge, subset = delta <= 2) cat.all.model.avg summary(cat.all.model.avg) #results are the same because symbiont density remains insig. ##pox is next## pox=data[,c(1,4,8:10,15)] ##test for outliers and normality## rosnerTest(pox$pox_stand) ##one pox <- pox[!(row.names(pox) %in% c("RI55")),] shapiro.test(pox$pox_stand) ##normal## ##first model with symbiont density continuous## pox_continuous=pox[,c(2:6)] pox_continuous = na.omit(pox_continuous) ##set up for dredging pox_global_model<- glm(pox_stand ~ . + .^4, data =pox_continuous, na.action = "na.fail") #dredge all models pox.all.dredge <- dredge(pox_global_model, beta = c("none"), evaluate = TRUE, rank = "AICc", cluster = "NULL") #view poxalase model output pox.all.dredge summary(pox.all.dredge) #average models between 0-2 delta AICc pox.all.model.avg <- model.avg(pox.all.dredge, subset = delta <= 2) pox.all.model.avg summary(pox.all.model.avg) #Model-averaged coefficients: #(full average) #Estimate Std. Error Adjusted SE z value Pr(>|z|) #(Intercept) 1.962535 0.383062 0.403102 4.869 1.1e-06 *** #symbiont_stand -0.001447 0.001581 0.001635 0.885 0.374 #(conditional average) #Estimate Std. Error Adjusted SE z value Pr(>|z|) #(Intercept) 1.962535 0.383062 0.403102 4.869 1.12e-06 *** #symbiont_stand -0.002476 0.001316 0.001423 1.739 0.0812 . #--- #Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 ##then model with poxegorical symbiont density## pox_bi=pox[,c(1:4,6)] pox_bi = na.omit(pox_bi) ##set up for dredging pox_global_model<- glm(pox_stand ~ . + .^4, data =pox_bi, na.action = "na.fail") #dredge all models pox.all.dredge <- dredge(pox_global_model, beta = c("none"), evaluate = TRUE, rank = "AICc", cluster = "NULL") #view poxalase model output pox.all.dredge summary(pox.all.dredge) #average models between 0-2 delta AICc pox.all.model.avg <- model.avg(pox.all.dredge, subset = delta <= 2) pox.all.model.avg summary(pox.all.model.avg) #Model-averaged coefficients: #(full average) #Estimate Std. Error Adjusted SE z value Pr(>|z|) #(Intercept) 1.2528 0.4788 0.5007 2.502 0.0123 * #typeWhite 0.6579 0.6223 0.6435 1.022 0.3066 #(conditional average) #Estimate Std. Error Adjusted SE z value Pr(>|z|) #(Intercept) 1.2528 0.4788 0.5007 2.502 0.0123 * #typeWhite 1.0062 0.4916 0.5318 1.892 0.0585 . #--- #Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 ##symbiont type near sig## ##dt is next## dt=data[,c(1,6,8:10,15)] ##test for outliers and normality## rosnerTest(dt$dt) ##none shapiro.test(dt$dt) ##normal## ##first model with symbiont density continuous## dt_continuous=dt[,c(2:6)] dt_continuous = na.omit(dt_continuous) ##set up for dredging dt_global_model<- glm(dt ~ . + .^4, data =dt_continuous, na.action = "na.fail") #dredge all models dt.all.dredge <- dredge(dt_global_model, beta = c("none"), evaluate = TRUE, rank = "AICc", cluster = "NULL") #view dtalase model output dt.all.dredge summary(dt.all.dredge) #average models between 0-2 delta AICc dt.all.model.avg <- model.avg(dt.all.dredge, subset = delta <= 2) dt.all.model.avg summary(dt.all.model.avg) #Model-averaged coefficients: #(full average) #Estimate Std. Error Adjusted SE z value Pr(>|z|) #(Intercept) 3.3547 0.1866 0.1955 17.156 <2e-16 *** #carb_stand 13.5320 20.3281 20.9210 0.647 0.518 #lipid_stand 5.3217 14.9541 15.5107 0.343 0.732 #(conditional average) #Estimate Std. Error Adjusted SE z value Pr(>|z|) #(Intercept) 3.3547 0.1866 0.1955 17.156 <2e-16 *** #carb_stand 32.6855 19.2896 20.7644 1.574 0.115 #lipid_stand 27.2706 23.3963 25.1851 1.083 0.279 #--- #Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 ##then model with categorical symbiont density## dt_bi=dt[,c(1:4,6)] dt_bi = na.omit(dt_bi) ##set up for dredging dt_global_model<- glm(dt ~ . + .^4, data =dt_bi, na.action = "na.fail") #dredge all models dt.all.dredge <- dredge(dt_global_model, beta = c("none"), evaluate = TRUE, rank = "AICc", cluster = "NULL") #view dtalase model output dt.all.dredge summary(dt.all.dredge) #average models between 0-2 delta AICc dt.all.model.avg <- model.avg(dt.all.dredge, subset = delta <= 2) dt.all.model.avg summary(dt.all.model.avg) #same because symbiont density isn't sig ##mel is next## mel=data[,c(1,7,8:10,15)] ##test for outliers and normality## rosnerTest(mel$mel_stand) ##one mel <- mel[!(row.names(mel) %in% c("RI51")),] shapiro.test(mel$mel_stand) ##normal## ##first model with symbiont density continuous## mel_continuous=mel[,c(2:6)] mel_continuous = na.omit(mel_continuous) ##set up for dredging mel_global_model<- glm(mel_stand ~ . + .^4, data =mel_continuous, na.action = "na.fail") #dredge all models mel.all.dredge <- dredge(mel_global_model, beta = c("none"), evaluate = TRUE, rank = "AICc", cluster = "NULL") #view melalase model output mel.all.dredge summary(mel.all.dredge) #only one sig model (carb +sym) summary(glm(mel_stand ~ carb_stand + symbiont_stand, data=mel_continuous)) #Coefficients: #Estimate Std. Error t value Pr(>|t|) #(Intercept) 2.609e-05 1.610e-05 1.621 0.12595 #carb_stand 5.439e-03 1.823e-03 2.984 0.00927 ** #symbiont_stand 1.966e-10 7.374e-11 2.667 0.01759 * ##symbionts and carbs are significant ##then model with categorical symbiont density## mel_bi=mel[,c(1:4,6)] mel_bi = na.omit(mel_bi) ##set up for dredging mel_global_model<- glm(mel_stand ~ . + .^4, data =mel_bi, na.action = "na.fail") #dredge all models mel.all.dredge <- dredge(mel_global_model, beta = c("none"), evaluate = TRUE, rank = "AICc", cluster = "NULL") #view melalase model output mel.all.dredge summary(mel.all.dredge) #average models between 0-2 delta AICc mel.all.model.avg <- model.avg(mel.all.dredge, subset = delta <= 2) mel.all.model.avg summary(mel.all.model.avg) #Model-averaged coefficients: #(full average) #Estimate Std. Error Adjusted SE z value Pr(>|z|) #(Intercept) 5.010e-05 1.915e-05 2.056e-05 2.436 0.0148 * #carb_stand 5.083e-03 2.114e-03 2.288e-03 2.222 0.0263 * #typeWhite -1.238e-05 1.769e-05 1.833e-05 0.675 0.4996 #(conditional average) #Estimate Std. Error Adjusted SE z value Pr(>|z|) #(Intercept) 5.010e-05 1.915e-05 2.056e-05 2.436 0.0148 * #carb_stand 5.083e-03 2.114e-03 2.288e-03 2.222 0.0263 * #typeWhite -2.746e-05 1.674e-05 1.821e-05 1.508 0.1316 #--- #Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 ##plot mel carb and mel symbiont stuff## ##fig1## melcarb=ggplot(mel, aes(mel_stand,carb_stand)) + geom_smooth(method = "lm", color = "black", lty = 1) + geom_point(aes(fill = type), pch = 21, color = "black", size = 3) + scale_fill_manual(values = c('Brown' = "burlywood4", 'White' = "white")) + theme_bw() + theme(panel.grid = element_blank()) + xlab("Melanin Concentration\n(mg melanin/mg tissue)") + ylab("Carbohydrate Concentration\n(mg lipid/mg tissue)") melsym=ggplot(mel, aes(mel_stand,symbiont_stand)) + geom_smooth(method = "lm", color = "black", lty = 1) + geom_point(aes(fill = type), pch = 21, color = "black", size = 3) + scale_fill_manual(values = c('Brown' = "burlywood4", 'White' = "white")) + theme_bw() + theme(panel.grid = element_blank()) + xlab("Melanin Concentration\n(mg melanin/mg tissue)") + ylab("Symbiont Density\n(cells/area)") ggarrange(melcarb, melsym, nrow = 1, ncol =2, align = "hv", common.legend = TRUE, legend="right") ##finally ppo ppo=data[,c(1,5,8:10,15)] ##test for outliers and normality## rosnerTest(ppo$ppo_stand) ##none shapiro.test(ppo$ppo_stand) ##non-normal## ppo$ppo_stand.norm <- log(ppo$ppo_stand) ##first model with symbiont density continuous## ppo_continuous=ppo[,c(3:7)] ppo_continuous = na.omit(ppo_continuous) ##set up for dredging ppo_global_model<- glm(ppo_stand.norm ~ . + .^4, data =ppo_continuous, na.action = "na.fail") #dredge all models ppo.all.dredge <- dredge(ppo_global_model, beta = c("none"), evaluate = TRUE, rank = "AICc", cluster = "NULL") #view ppoalase model output ppo.all.dredge summary(ppo.all.dredge) #average models between 0-2 delta AICc ppo.all.model.avg <- model.avg(ppo.all.dredge, subset = delta <= 2) ppo.all.model.avg summary(ppo.all.model.avg) #Model-averaged coefficients: #(full average) #Estimate Std. Error Adjusted SE z value Pr(>|z|) #(Intercept) -1.39411 0.14375 0.15292 9.117 <2e-16 *** #spermpolypstand 0.01346 0.02779 0.02893 0.465 0.642 #(conditional average) #Estimate Std. Error Adjusted SE z value Pr(>|z|) #(Intercept) -1.39411 0.14375 0.15292 9.117 <2e-16 *** #spermpolypstand 0.04062 0.03504 0.03772 1.077 0.282 #--- #Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 ##nothing sig ##then model with categorical symbiont density## ppo_bi=ppo[,c(1,3:5,7)] ppo_bi = na.omit(ppo_bi) ##set up for dredging ppo_global_model<- glm(ppo_stand.norm ~ . + .^4, data =ppo_bi, na.action = "na.fail") #dredge all models ppo.all.dredge <- dredge(ppo_global_model, beta = c("none"), evaluate = TRUE, rank = "AICc", cluster = "NULL") #view ppoalase model output ppo.all.dredge summary(ppo.all.dredge) #only one sig model, and its the null one.