R Notebook

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Cmd+Shift+Enter.

IsosWing1 <- read_csv("IsogenicWingArea_60LINES_omitted10lines.csv", 
     col_types = cols(Diet_dilution = col_double(), 
       Line = col_character(), Replicate = col_character(), Block = col_character(),
        Temp = col_character()))
head (IsosWing1)

nrow (IsosWing1)

[1] 3899

function to rotate text on the x and y axes

rotatedAxisElementText = function(angle,position='x'){
  angle     = angle[1]; 
  position  = position[1]
  positions = list(x=0,y=90,top=180,right=270)
  if(!position %in% names(positions))
    stop(sprintf("'position' must be one of [%s]",paste(names(positions),collapse=", ")),call.=FALSE)
  if(!is.numeric(angle))
    stop("'angle' must be numeric",call.=FALSE)
  rads  = (angle - positions[[ position ]])*pi/180
  hjust = 0.5*(1 - sin(rads))
  vjust = 0.5*(1 + cos(rads))
  element_text(angle=angle,vjust=vjust,hjust=hjust)
}

rr means1 <- IsosWing1 %>% group_by(Line, Temp, Diet_dilution) %>% summarise(mean = mean(wing_area)) means1

       mean
1 0.9371405

rr #write.csv(means1, _means.csv) #for orginial 60lines dataset new_dataframe_means <- IsosWing %>% group_by(Line, Temp, Diet_dilution) %>% mutate (means= mean(wing_area)) new_dataframe_means

# A tibble: 3,010 x 7
# Groups:   Line, Temp, Diet_dilution [245]
   Diet        Temp  Line Replicate wing_area Diet_dilution means
   <chr>      <int> <int>     <int>     <dbl>         <dbl> <dbl>
 1 2:3(X1)       25    29         1     1.05            100 0.930
 2 2:3(X1)       25    29         2     0.971           100 0.930
 3 2:3(X1)       25    29         3     1.03            100 0.930
 4 2:3(X1)       25    29         4     0.985           100 0.930
 5 2:3(X1)       25    29         5     0.985           100 0.930
 6 2:3(X1)       25    29         6     1.00            100 0.930
 7 2:3(X1)       25    29         7     0.99            100 0.930
 8 2:3(X1)       25    29         8     0.892           100 0.930
 9 2:3(X1)       25    29         9     0.983           100 0.930
10 2:3(X0.25)    25    29         1     0.946            25 0.930
# … with 3,000 more rows

rr means2 <- IsosWing %>% group_by(Line) %>% summarise(mean = mean(wing_area), sd = sd(wing_area)) means2

       mean       sd
1 0.9299233 0.113709

lrtest (IsosWing.lm4,IsosWing.lm2)

Likelihood ratio test

Model 1: wing_area ~ 1 + Diet_dilution * Temp + (1 | Line) + (1 | Diet_dilution:Line) + 
    (1 | Temp:Line) + (1 | Diet_dilution:Temp:Line) + (1 | Block)
Model 2: wing_area ~ -1 + Diet_dilution * Temp + (1 | Block) + (1 | Line)
  #Df LogLik Df  Chisq Pr(>Chisq)    
1  10 4727.0                         
2   7 4484.4 -3 485.07  < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

chi2 <- 2*(summary(mod5)$logLik - summary(mod4)$logLik)
1-pchisq(chi2, 3)

'log Lik.' 3.026083e-10 (df=16)

rr summary(mod2b)

#making a new diet_dilution by temp to see if the model fits better IsosWing\(diet_dilution_by_temp<-as.factor(paste(IsosWing\)Diet_dilution,IsosWing$Temp,sep=;)) mod3<-lme(wing_area ~ 1+Diet_dilution * Temp,random=c(~-1+diet_dilution_by_temp|Line),data=IsosWing,control=list(opt=)) summary (mod3)

anova(mod3,mod2b) anova(mod2b,mod2a)

summary(mod3) confint(mod3) ranef(model3.1)

library() library(4) dotplot(ranef(mod3, condVar=TRUE), lattice.options=list(layout=c(1,2))) mod3<-lmer(wing_area ~ 1+Diet_dilution * Temp + (-1+diet_dilution|Line), data=IsosWing) dotplot(ranef(mod3, condVar=TRUE), lattice.options=list(layout=c(1,2)))

IsosWing_scaled <- IsosWing1 %>% mutate(Diet_dilution_scaled = scale(Diet_dilution, center = T, scale = T), Diet_Temp = paste(Diet_dilution, Temp, sep = _)) IsosWing_scaled

model3.1 <- lmer(wing_area ~ 1+Diet_dilution_scaled * Temp + (1+Diet_dilution_scaled*Temp|Line) + (1|Block), REML = FALSE, data=IsosWing_scaled) plot(model3.1)

summary(model3.1) car::Anova(model3.1) model3.2 <- lmer(wing_area ~ 1+Diet_dilution_scaled + Temp + (1+Diet_dilution_scaled*Temp|Line) + (1|Block), REML = FALSE, data=IsosWing_scaled)

summary (model3.2) dotplot(ranef(model3.1, condVar=TRUE), lattice.options=list(layout=c(1,2)))

#ranef gives us the BLUPs, which represents the response of a given genotype to the fixed effect of the diet at each temp, as difference betwen the genotype’s predicted response and the population level average predicted response genotype_blups <- ranef(model3.1)$Line

genotype_index <- as.factor(genotype_blups$Line) genotype_index <- data.frame(rownames(genotype_blups), genotype_blups) genotype_data <- cbind( genotype_index, genotype_blups) colnames(genotype_data) <- c(, _int, _slope, ,_slope_by_temp )

#there are two random effects for RRMs: intercept and slopes #BLUPs by intercept, indicates the difference in genotype elevation relative to the population average, more +ve BLUP intercept means that the genotype’s reaction norm occurs above the pop-level average response. -ve indicated the genotype’s reaction norm is below pop-level average response. these are not a measure of plasticity. ggplot(genotype_data, aes(Line, BLUP_int )) + geom_point(aes(group = Line, colour = Line), size = 4) + ylab(intercept estimate) + geom_hline(yintercept = 0, lty = 2) + theme_classic()

#BLUPs by slope estimate is the difference in slope(relative steepness of change) between the pop-level average response and the response of the individual genotype. ggplot(genotype_data, aes(Line, BLUP_slope)) + geom_point(aes(group = Line, colour = Line), size = 4) + ylab((BLUP slope estimate)) + geom_hline(yintercept = 0, lty = 2) + theme_classic()

#is there any coorelation with BLUP intercept and BLUP slope? ggplot(genotype_data, aes(BLUP_int, BLUP_slope)) + geom_point(aes(group = Line, colour = Line), size = 4) + xlab(intercept estimate) + ylab(slope estimate) + theme_classic() #answer: not really a positive corelation, as can be seen from summary table of the model as well?

#rank BLUPs in order, sort by slope of most to least plastic. more negative BLUP slope estimate represent greater plasticity and more positive BLUP slope estimates represents less plasticity in the relative wing area across nutritional conditions at 25 vs 28. genotype_data\(genotype_ordered <- factor(genotype_data\)Line, levels = genotype_data\(Line[order(genotype_data\)BLUP_slope_by_temp)]) plotf<-ggplot(genotype_data, aes(genotype_ordered, BLUP_slope_by_temp)) + geom_bar(stat = , aes(group = Line, fill = Line)) + xlab((in ranked order of plasticity)) + ylab((BLUP slope estimate)) + theme_classic()

pdf(./plotf.pdf, useDingbats=FALSE, width = 12 , height = 7 ) grid.arrange(plotf) dev.off()

rr #plot the raw values of wing area against mean-centred and scaled diet_dilution for each line. ggplot(IsosWing_scaled, aes(x = Diet_dilution_scaled, y = wing_area, group = Line)) + geom_line(aes(colour = Line)) + ylab(wing area) + xlab(-centered and scaled rearing diet_dilution) + facet_grid(~ Temp)+theme_classic()

#fit the model on this to visualise fit?

pdf("./50lines_population_mean_response.pdf", useDingbats=FALSE) 
grid.arrange(PLOT)
dev.off()

null device 
          1 

pdf("./50lines_response.pdf", useDingbats=FALSE) 
grid.arrange(PLOT)
dev.off()

IsosWing1$Temp <- as.numeric(IsosWing1$Temp)
IsosWing1$Diet_dilution <- as.numeric(IsosWing1$Diet_dilution)
as.numeric(as.character(Temp))
IsosDT.lm <- lmer (wing_area ~ Temp * Line * Diet_dilution + (1|Block), data = IsosWing1) #line as fixed effect since we have chosen them from prior analysis (RRMs)
anova (IsosDT.lm) #supplimentary table with wing area and dev time, but refer to it saying evrything differs
Isos_trend1 <- emtrends(IsosDT.lm, ~ Temp | Line, var = "Diet_dilution")
multcomp::cld(Isos_trend1)

Isos_trend2 <- emtrends(IsosDT.lm, ~ Line | Temp, var = "Diet_dilution")
multcomp::cld(Isos_trend2)

Isos_trend3 <- emtrends(IsosDT.lm, ~ Diet_dilution | Line, var = "Temp")
multcomp::cld(Isos_trend3)

New$Temp <- as.numeric(New$Temp)
New$Diet_dilution <- as.numeric(New$Diet_dilution)
as.numeric(as.character(Temp))
IsosDT.lm <- lmer (wing_area ~ Temp * Line * Diet_dilution + (1|Block), data = New) #line as fixed effect since we have chosen them from prior analysis (RRMs)
anova (IsosDT.lm) #supplimentary table with wing area and dev time, but refer to it saying evrything differs
Isos_trend1 <- emtrends(IsosDT.lm, ~ Temp | Line, var = "Diet_dilution")
multcomp::cld(Isos_trend1)

Isos_trend2 <- emtrends(IsosDT.lm, ~ Line | Temp, var = "Diet_dilution")
multcomp::cld(Isos_trend2)

Isos_trend3 <- emtrends(IsosDT.lm, ~ Diet_dilution | Line, var = "Temp")
multcomp::cld(Isos_trend3)

Isos_trend1 <- emtrends(IsosDT.lm, ~ Temp | Line, var = "Diet_dilution_scaled")

Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'pbkrtest.limit = 3899' (or larger)
[or, globally, 'set emm_options(pbkrtest.limit = 3899)' or larger];
but be warned that this may result in large computation time and memory use.
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'lmerTest.limit = 3899' (or larger)
[or, globally, 'set emm_options(lmerTest.limit = 3899)' or larger];
but be warned that this may result in large computation time and memory use.

rr pdf(./5lines_response.pdf, useDingbats=FALSE) grid.arrange(PLOT2) dev.off()

null device 
          1 

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Cmd+Option+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Cmd+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.