#### Costs of development and temperature (meta-analysis 2) #### library(asreml) library(ape) library(MCMCglmm) library(ggplot2) x<-rnorm(100) y<-rnorm(100) m1<-asreml(y~x) pin <- dget("pin.R") #put into same location as asreml # Data data <- read.csv("matched data2.csv") data$Species <- as.factor(data$Species) # Tree tree <- read.tree("matched tree2.phy") tree$node.label <- c(1:length(tree$node.label)) par(mfrow = c(1,1), mar = rep(0,4)) plot(tree, cex = 0.5) par(mfrow = c(2,2)) plot(propn_cost_warm_dividedby_mid ~ deta_T_range, data = data) boxplot(data$propn_cost_warm_dividedby_mid) # Convert the tree into the matrix format expected by ASReml data$animal <- as.factor(data$Species) data.IA <- inverseA(tree, nodes = "TIPS") data.IAasreml <- sm2asreml(data.IA$Ainv, data.IA$node.names) data$animal <- as.character(data$animal) # Check for missing data missing.from.data <- setdiff(tree$tip.label, data$animal) missing.from.data missing.from.tree <- setdiff(data$animal, tree$tip.label) missing.from.tree # Run a model m1 <- asreml(fixed = propn_cost_warm_dividedby_mid ~ deta_T_range, random= ~giv(animal,var=T), data = data, ginverse = list(animal=data.IAasreml)) summary(m1) # $call # asreml(fixed = propn_cost_warm_dividedby_mid ~ deta_T_range, # random = ~giv(animal, var = T), data = data, ginverse = list(animal = data.IAasreml)) # # $loglik # [1] 81.9808 # # $nedf # [1] 70 # # $sigma # [1] 0.1741034 # # $varcomp # gamma component std.error z.ratio constraint # giv(animal, var = T).giv 0.2502187 0.007584625 0.009122235 0.8314437 Positive # R!variance 1.0000000 0.030311987 0.005749083 5.2724911 Positive # # attr(,"class") # [1] "summary.asreml" summary(m1,all=T)$coef.fi #PMM fixed effects # solution std error z ratio # deta_T_range 0.1407479 0.19281954 0.7299461 # (Intercept) 0.7626932 0.06367108 11.9786450 wald.asreml(m1, ssType="conditional", denDF="numeric") #Significance of fixed effects # LogLik S2 DF wall cpu # 81.9808 0.0303 70 08:48:48 0.0 # 81.9808 0.0303 70 08:48:48 0.0 # 81.9808 0.0303 70 08:48:48 0.0 # 81.9808 0.0303 70 08:48:48 0.0 # # LogLikelihood Converged # $Wald # Df denDF F.inc F.con Margin Pr # (Intercept) 1 4.0 217.3000 143.3000 0.0002941669 # deta_T_range 1 69.8 0.5337 0.5337 A 0.4675086334 # # $stratumVariances # NULL pin(m1, ~ V1/(V1+V2)) #Calculate phylogenetic heritability # Estimate SE # 0.2001399 0.2048172 m1.lm <- lm(propn_cost_warm_dividedby_mid ~ deta_T_range, data = droplevels(subset(data, !is.na(deta_T_range) & !is.na(animal)))) summary(m1.lm) # Call: # lm(formula = propn_cost_warm_dividedby_mid ~ deta_T_range, data = droplevels(subset(data, # !is.na(deta_T_range) & !is.na(animal)))) # # Residuals: # Min 1Q Median 3Q Max # -0.57255 -0.12037 0.05217 0.13760 0.31884 # # Coefficients: # Estimate Std. Error t value Pr(>|t|) # (Intercept) 0.78822 0.04105 19.199 <2e-16 *** # deta_T_range 0.01619 0.19372 0.084 0.934 # --- # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 # # Residual standard error: 0.1849 on 70 degrees of freedom # Multiple R-squared: 9.974e-05, Adjusted R-squared: -0.01418 # F-statistic: 0.006982 on 1 and 70 DF, p-value: 0.9336 m1.lm.predictions <- predict(m1.lm, newdata = data.frame(deta_T_range = seq(from = min(data$deta_T_range),to = max(data$deta_T_range),by = 0.1)),interval = "confidence") m1.lm.predictions #Run a model m1 <- asreml(fixed = propn_cost_warm_dividedby_mid ~ deta_T_range, random= ~giv(animal,var=T) + Species, data = data, ginverse = list(animal=data.IAasreml)) summary(m1) # $call # asreml(fixed = propn_cost_warm_dividedby_mid ~ deta_T_range, # random = ~giv(animal, var = T) + Species, data = data, ginverse = list(animal = data.IAasreml)) # # $loglik # [1] 81.79135 # # $nedf # [1] 70 # # $sigma # [1] 0.1706793 # # $varcomp # gamma component std.error z.ratio constraint # giv(animal, var = T).giv 0.100000000 0.0029131414 0.004560806 0.638733951 Singular # Species!Species.var 0.006324555 0.0001842432 0.051880284 0.003551315 ? # R!variance 1.000000000 0.0291314144 0.045608057 0.638733951 Positive # # attr(,"class") # [1] "summary.asreml" summary(m1,all=T)$coef.fi #PMM fixed effects # solution std error z ratio # deta_T_range 0.09200494 0.19330014 0.4759693 # (Intercept) 0.77374377 0.05100348 15.1704117 wald.asreml(m1, ssType="conditional", denDF="numeric") #Significance of fixed effects # LogLik S2 DF wall cpu # 81.8189 0.0317 70 08:52:35 0.0 (1 restrained) # 1 singularities inAverage Information Matrix # # $Wald # Df denDF F.inc F.con Margin Pr # (Intercept) 1 0.9 421.8000 223.9000 0.0583671 # deta_T_range 1 0.9 0.2468 0.2481 A 0.7148210 # # $stratumVariances # NULL pin(m1, ~ V1/(V1+V2+V3)) # Estimate SE # 0.09038939 0.1426045 #Plot the data par(mfrow = c(1,1), mar = c(4,4,0,0)+0.5) plot(propn_cost_warm_dividedby_mid ~ deta_T_range, data = data) m2 <- asreml(fixed = propn_cost_warm_dividedby_mid ~ deta_T_range, random= ~giv(animal, var=T) + Species, data = data, #weights = "sqrt.n", ginverse = list(animal=data.IAasreml)) summary(m2) summary(m2,all=T)$coef.fi #PMM fixed effects wald.asreml(m2, ssType="conditional", denDF="numeric") #Significance of fixed effects pin(m2, ~ V1/(V1+V2+V3)) data2 <- subset(data, !is.na(animal)) m2.without.phylo <- asreml(fixed = propn_cost_warm_dividedby_mid ~ deta_T_range, random= ~Species, #weights = "sqrt.n" data = data) summary(m2.without.phylo) #Significance of random effects of phylogeny 2*(m2.without.phylo$loglik-m2$loglik) #ChiSq 1-pchisq(2*(m2.without.phylo$loglik-m2$loglik),1) #p #Plot the data par(mfrow = c(1,1), mar = c(4,4,0,0)+0.5) plot(propn_cost_warm_dividedby_mid ~ deta_T_range, ylim = c(0.1, 1.2), data = data) #type = "n") # plot(Delta.Cost ~ Delta.T.Range, data = data.2) polygon(x = c(seq(from = min(data$deta_T_range), to = max(data$deta_T_range), by = 0.1), rev(seq(from = min(data$deta_T_range), to = max(data$deta_T_range), by = 0.1))), y = c(m1.lm.predictions[,2], rev(m1.lm.predictions[,3])), col = "grey80", border = NA) lines(x = seq(from = min(data$deta_T_range), to = max(data$deta_T_range), by = 0.1), y = m1.lm.predictions[,1]) abline(h = 1, lty = 3) points(propn_cost_warm_dividedby_mid ~ deta_T_range, data = data, pch = 21, col = "white", bg = "black") # Boxplot figure (Figure 4) data <- read.csv("matched data2.csv") boxplot(data$propn_cost_warm_dividedby_mid) + abline(h = 1, lty = 3) p <- ggplot(data, aes(x = group, y = propn_cost_warm_dividedby_mid, fill = group)) + geom_boxplot(alpha = 0.80) + geom_point(aes(fill = group), alpha = 0.3, size = 3, shape = 21, position = position_jitterdodge()) + scale_y_continuous(limits=c(0,1.3), breaks=seq(0.2,1.2,0.4)) + theme_bw() + abline(h=1, lwd = 3, col = 'black') + theme(text = element_text(size = 10), axis.title.x = element_blank(), axis.title.y = element_blank(), panel.grid = element_blank(), legend.position = "none") p + geom_hline(aes(yintercept=1, linetype="dashed"))