# Evolution of sexual conflict - Single mating # Analysis of male harm: Effects of focal males on tester females' fitness. # Experiment carried out at generation 32. # Response variables: LRS (Lifetime Reproductive Success: number of adults produced by each female during her whole life), # longev: longevity, in days. # Factors: treatment, ss, structure, gen (generation), # Random factor: line (unique code 1 to 16). # maleid: Unique code for each indididual # Covariables: fel and mel: female and male elytrum length (proxy for body size). Always measured by Miguel Lozano. #Load libraries library(lme4) library(car) library(MASS) # Input data data = read.table("focalmalessinglecop.txt",header=TRUE,sep="\t") dim(data) names(data) head(data) str(data) # Define factors data$treatment <- as.factor(data$treatment) data$ss <- as.factor(data$ss) data$structure <- as.factor(data$structure) data$str <- as.factor(data$str) data$line <- as.factor(data$line) data$maleid<- as.factor(data$maleid) data$malenumber<- as.factor(data$malenumber) str(data) #Balanced dataset: table(data$str, data$ss) #Standard code for basic statistics and initial inspection of distributions, etc. not included. # Centering covariates data$cent.fel<-(data$fel-mean(data$fel,na.rm=TRUE)) hist(data$cent.fel) mean(data$cent.fel, na.rm=T) shapiro.test(data$cent.fel) data$cent.mel<-(data$mel-mean(data$mel,na.rm=TRUE)) hist(data$cent.mel) mean(data$cent.mel, na.rm=T) shapiro.test(data$cent.mel) data$cent.longev<-(data$longev-mean(data$longev,na.rm=TRUE)) hist(data$cent.longev) mean(data$cent.longev, na.rm=T) shapiro.test(data$cent.longev) data$cent.LRS<-(data$LRS-mean(data$LRS,na.rm=TRUE)) hist(data$cent.LRS) mean(data$cent.LRS, na.rm=T) shapiro.test(data$cent.LRS) #### LONGEVITY LMM with Random intercept and Random slope#### # Model estimating random intercept, random slopes, and the covariance between slope and intercept. M1<-lmer(longev~ss+str+cent.fel+cent.mel+cent.LRS+ ss:str+ ss:cent.fel+ss:cent.mel+ss:cent.LRS+ str:cent.fel+str:cent.mel+str:cent.LRS+ (1+cent.fel+cent.mel+cent.LRS|line), REML=FALSE, data = data) summary(M1) Anova(M1) #Chi square. Wald test. Type II # This model has got problems of convergence. Barr et al. (2013.Journal of Memory and Language 68:255-278) # conclude the following "(the results of their simulations)...This suggests that when maximal LMEMs fail to converge, # dropping within-unit random intercepts or random correlations are both viable options for simplifying the # random effects structure”. # In other words, apparently the decision to model the correlations in the random structure doesn’t seem to inflate type I error. # So, below I run a model without the correlation structure, and it doesn't present convergence issues: M1<-lmer(longev~ss+str+cent.fel+cent.mel+cent.LRS+ ss:str+ ss:cent.fel+ss:cent.mel+ss:cent.LRS+ str:cent.fel+str:cent.mel+str:cent.LRS+ (1|line)+(0+cent.fel|line)+(0+cent.mel|line)+(0+cent.LRS|line), REML=FALSE, data = data) summary(M1) Anova(M1) #Chi square. Wald test. Type II #Final model with REML M1<-lmer(longev~ss+str+cent.fel+cent.mel+cent.LRS+ ss:str+ ss:cent.fel+ss:cent.mel+ss:cent.LRS+ str:cent.fel+str:cent.mel+str:cent.LRS+ (1|line)+(0+cent.fel|line)+(0+cent.mel|line)+(0+cent.LRS|line), REML=TRUE, data = data) summary(M1) # Diagnostic plots: par(mfrow = c(2, 2)) plot(M1) diagnostics.plot(M1) #Checking basic assumptions #### LRS models #### #### LRS LMM with Random intercept and Random slope #### #Based on diagnostic plots LRS may need transformation for better fit par(mfrow = c(1, 1)) hist(data$LRS) data$sqLRS <- ((data$LRS)^2) # Squared hist(data$sqLRS) #This is the one providing better fit to normal curve. # Model M1<-lmer(sqLRS~ss+str+cent.fel+cent.mel+cent.longev+ ss:str+ ss:cent.fel+ss:cent.mel+ss:cent.longev+ str:cent.fel+str:cent.mel+str:cent.longev+ (1|line)+(0+cent.fel|line)+(0+cent.mel|line)+(0+cent.longev|line), REML=FALSE, data = data) summary(M1) Anova(M1) #Chi square. Wald test. Type II #Final model with REML M1<-lmer(sqLRS~ss+str+cent.fel+cent.mel+cent.longev+ ss:str+ ss:cent.fel+ss:cent.mel+ss:cent.longev+ str:cent.fel+str:cent.mel+str:cent.longev+ (1|line)+(0+cent.fel|line)+(0+cent.mel|line)+(0+cent.longev|line), REML=TRUE, data = data) summary(M1) # Diagnostic plots: par(mfrow = c(2, 2)) plot(M1) diagnostics.plot(M1) #Checking basic assumptions