# WT vs IRF3IRF7 (Fig. 1)
library(lme4)
## Loading required package: Matrix
library(ggplot2)
library(readxl)
library(readr)
library(reshape2)
library(tidyr)
##
## Attaching package: 'tidyr'
## The following object is masked from 'package:reshape2':
##
## smiths
## The following objects are masked from 'package:Matrix':
##
## expand, pack, unpack
library(plyr)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:plyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
d <- read_excel("C:/Users/kmaciag/Desktop/IRF3 Paper - Local/2013-05-03 WT vs IRF3-IRF7 DKO anti-IFNAR pre or post and RNA-seq - full data.xlsx", sheet = "DKO vs WT - for GCA")
# Remove un-needed timepoints to standardize
d$Timept <- as.numeric(gsub("[^0-9.]", "", d$Timept))
d <- d[d$Timept!=60,]
# would prob keep 6, 14, 26...
# discard 20, ?32, 50
# melt data. no normalizing to BL
d2 <- d %>% gather(TechRep, lumi, -Timept, -Genotype, -IFNgTiming, -IFNgStim)
# log transform
d2$log10bact <- log(d2$lumi,10) # set Y value as log2 luminescence
# plot log-transformed
ggplot(d2, aes(Timept, log10bact, color=Genotype)) + geom_point()+
theme_minimal(base_size=20) +
facet_grid(IFNgTiming ~ IFNgStim)
#########################
# Cycle through the facets
# d2_facet <-d2[d2$IFNgTiming=="prestimIFNg" & d2$IFNgStim=="0U",]
# d2_facet <-d2[d2$IFNgTiming=="prestimIFNg" & d2$IFNgStim=="3U",]
# d2_facet <-d2[d2$IFNgTiming=="prestimIFNg" & d2$IFNgStim=="10U",]
# d2_facet <-d2[d2$IFNgTiming=="prestimIFNg" & d2$IFNgStim=="100U",]
d2_facet <-d2[d2$IFNgTiming=="postIFNg" & d2$IFNgStim=="0U",]
# d2_facet <-d2[d2$IFNgTiming=="postIFNg" & d2$IFNgStim=="3U",]
# d2_facet <-d2[d2$IFNgTiming=="postIFNg" & d2$IFNgStim=="10U",]
# d2_facet <-d2[d2$IFNgTiming=="postIFNg" & d2$IFNgStim=="100U",]
##### GCA - Quadratic ######
d3 <- d2_facet
# add timeBin
orderedTimepts <- unique(d3$Timept)[order(unique(d3$Timept))]
d3$timeBin <- match(d3$Timept,orderedTimepts)
# add 3rd degree orthogonal polynomial
t <- poly(unique(d3$Timept), 2) # one row for each timepoint, two cols for three dimensions
d3[,paste("ot", 1:2, sep="")] <- t[d3$timeBin, 1:2] # placeholders in ot1, ot2 for now
# Create models
# no group effects model
# poly_m.base <- lmer(log10bact ~ (ot1+ot2) + (1+ot1+ot2|BiolRepl/TechRep), data=d3, REML = FALSE)
poly_m.base <- lmer(log10bact ~ (ot1+ot2) + (1+ot1+ot2|TechRep), data=d3, REML = FALSE)
## boundary (singular) fit: see ?isSingular
# add effect of Genotype on intercept
# poly_m.0 <- lmer(log10bact ~ (ot1+ot2) + Genotype + (1+ot1+ot2|BiolRepl/TechRep), data=d3, REML = FALSE)
poly_m.0 <- lmer(log10bact ~ (ot1+ot2) + Genotype + (1+ot1+ot2|TechRep), data=d3, REML = FALSE)
## boundary (singular) fit: see ?isSingular
# add effect of Genotype on linear term
# poly_m.1 <- lmer(log10bact ~ (ot1+ot2) + Genotype + ot1:Genotype + (1+ot1+ot2|BiolRepl/TechRep), data=d3, REML = FALSE)
poly_m.1 <- lmer(log10bact ~ (ot1+ot2) + Genotype + ot1:Genotype + (1+ot1+ot2|TechRep), data=d3, REML = FALSE)
## boundary (singular) fit: see ?isSingular
# add effect of Genotype on quadratic term
# poly_m.2 <- lmer(log10bact ~ (ot1+ot2) + Genotype + ot1:Genotype + ot2:Genotype + (1+ot1+ot2|BiolRepl/TechRep), data=d3, REML = FALSE)
poly_m.2 <- lmer(log10bact ~ (ot1+ot2) + Genotype + ot1:Genotype + ot2:Genotype + (1+ot1+ot2|TechRep), data=d3, REML = FALSE)
## boundary (singular) fit: see ?isSingular
# compare the goodness of fit of the models
anova(poly_m.base, poly_m.0, poly_m.1, poly_m.2)
## Data: d3
## Models:
## poly_m.base: log10bact ~ (ot1 + ot2) + (1 + ot1 + ot2 | TechRep)
## poly_m.0: log10bact ~ (ot1 + ot2) + Genotype + (1 + ot1 + ot2 | TechRep)
## poly_m.1: log10bact ~ (ot1 + ot2) + Genotype + ot1:Genotype + (1 + ot1 +
## poly_m.1: ot2 | TechRep)
## poly_m.2: log10bact ~ (ot1 + ot2) + Genotype + ot1:Genotype + ot2:Genotype +
## poly_m.2: (1 + ot1 + ot2 | TechRep)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## poly_m.base 10 108.66 130.25 -44.331 88.662
## poly_m.0 11 110.38 134.12 -44.187 88.375 0.2874 1 0.5919
## poly_m.1 12 112.36 138.26 -44.179 88.358 0.0169 1 0.8967
## poly_m.2 13 114.29 142.36 -44.145 88.291 0.0668 1 0.7961
#Add model fits to d3 dataframe
d3$base <- predict(poly_m.base)
d3$m0 <- predict(poly_m.0)
d3$m1 <- predict(poly_m.1)
d3$m2 <- predict(poly_m.2)
## Show models on plots
# show info about base model fit
summary(poly_m.base)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: log10bact ~ (ot1 + ot2) + (1 + ot1 + ot2 | TechRep)
## Data: d3
##
## AIC BIC logLik deviance df.resid
## 108.7 130.3 -44.3 88.7 54
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.4129 -0.4797 0.2600 0.6234 1.4286
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## TechRep (Intercept) 0.000e+00 0.000e+00
## ot1 6.780e-11 8.234e-06 NaN
## ot2 1.111e-11 3.333e-06 NaN 0.80
## Residual 2.340e-01 4.837e-01
## Number of obs: 64, groups: TechRep, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 5.11991 0.06046 84.677
## ot1 0.98206 0.17102 5.742
## ot2 0.52365 0.17102 3.062
##
## Correlation of Fixed Effects:
## (Intr) ot1
## ot1 0.000
## ot2 0.000 0.000
## convergence code: 0
## boundary (singular) fit: see ?isSingular
# get parameter estimates and estimate p-values
m.coefs <- data.frame(coef(summary(poly_m.base)))
m.coefs$p <- 2*(1-pnorm(abs(m.coefs$t.value)))
m.coefs
## Estimate Std..Error t.value p
## (Intercept) 5.1199099 0.06046394 84.677074 0.000000e+00
## ot1 0.9820562 0.17101786 5.742419 9.333341e-09
## ot2 0.5236483 0.17101786 3.061951 2.198996e-03
# Plot base model predictions
ggplot(d3, aes(x = Timept, y = base, color = Genotype)) + geom_point() + stat_summary(fun.data=mean_se, fun = mean, geom = 'line', aes(group = Genotype)) + theme_minimal(base_size = 20)
# Overlay base model over summary data
ggplot(d3, aes(Timept, log10bact, color=Genotype)) +
stat_summary(fun.data=mean_se, geom="pointrange") +
stat_summary(aes(y=base, linetype=Genotype), fun=mean, geom="line")+theme_minimal(base_size = 20)
# show info about m0 model fit
summary(poly_m.0)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: log10bact ~ (ot1 + ot2) + Genotype + (1 + ot1 + ot2 | TechRep)
## Data: d3
##
## AIC BIC logLik deviance df.resid
## 110.4 134.1 -44.2 88.4 53
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.4185 -0.5413 0.1935 0.6545 1.3647
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## TechRep (Intercept) 0.000e+00 0.000e+00
## ot1 6.700e-11 8.185e-06 NaN
## ot2 1.124e-11 3.352e-06 NaN 0.80
## Residual 2.329e-01 4.826e-01
## Number of obs: 64, groups: TechRep, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 5.15229 0.08532 60.390
## ot1 0.98206 0.17063 5.755
## ot2 0.52365 0.17063 3.069
## GenotypeWT -0.06476 0.12066 -0.537
##
## Correlation of Fixed Effects:
## (Intr) ot1 ot2
## ot1 0.000
## ot2 0.000 0.000
## GenotypeWT -0.707 0.000 0.000
## convergence code: 0
## boundary (singular) fit: see ?isSingular
# get parameter estimates and estimate p-values
m.coefs <- data.frame(coef(summary(poly_m.0)))
m.coefs$p <- 2*(1-pnorm(abs(m.coefs$t.value)))
m.coefs
## Estimate Std..Error t.value p
## (Intercept) 5.15229031 0.08531712 60.3898726 0.000000e+00
## ot1 0.98205625 0.17063425 5.7553290 8.647321e-09
## ot2 0.52364827 0.17063425 3.0688345 2.148956e-03
## GenotypeWT -0.06476083 0.12065664 -0.5367366 5.914496e-01
# Plot m0 model predictions
ggplot(d3, aes(x = Timept, y = m0, color = Genotype)) + geom_point() + stat_summary(fun.data=mean_se, fun = mean, geom = 'line', aes(group = Genotype)) + theme_minimal(base_size = 20)
# Overlay m0 model over summary data
ggplot(d3, aes(Timept, log10bact, color=Genotype)) +
stat_summary(fun.data=mean_se, geom="pointrange") +
stat_summary(aes(y=m0, linetype=Genotype), fun=mean, geom="line")+theme_minimal(base_size = 20)
# show info about m1 model fit
summary(poly_m.1)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: log10bact ~ (ot1 + ot2) + Genotype + ot1:Genotype + (1 + ot1 +
## ot2 | TechRep)
## Data: d3
##
## AIC BIC logLik deviance df.resid
## 112.4 138.3 -44.2 88.4 52
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.4019 -0.5400 0.1879 0.6523 1.3910
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## TechRep (Intercept) 0.000e+00 0.000e+00
## ot1 5.890e-13 7.674e-07 NaN
## ot2 8.083e-13 8.991e-07 NaN 1.00
## Residual 2.329e-01 4.826e-01
## Number of obs: 64, groups: TechRep, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 5.15229 0.08531 60.398
## ot1 1.00421 0.24128 4.162
## ot2 0.52365 0.17061 3.069
## GenotypeWT -0.06476 0.12064 -0.537
## ot1:GenotypeWT -0.04431 0.34122 -0.130
##
## Correlation of Fixed Effects:
## (Intr) ot1 ot2 GntyWT
## ot1 0.000
## ot2 0.000 0.000
## GenotypeWT -0.707 0.000 0.000
## ot1:GntypWT 0.000 -0.707 0.000 0.000
## convergence code: 0
## boundary (singular) fit: see ?isSingular
# get parameter estimates and estimate p-values
m.coefs <- data.frame(coef(summary(poly_m.1)))
m.coefs$p <- 2*(1-pnorm(abs(m.coefs$t.value)))
m.coefs
## Estimate Std..Error t.value p
## (Intercept) 5.15229031 0.08530589 60.3978284 0.000000e+00
## ot1 1.00421198 0.24128148 4.1619936 3.154813e-05
## ot2 0.52364827 0.17061177 3.0692387 2.146050e-03
## GenotypeWT -0.06476083 0.12064074 -0.5368073 5.914008e-01
## ot1:GenotypeWT -0.04431146 0.34122355 -0.1298605 8.966768e-01
# Plot m1 model predictions
ggplot(d3, aes(x = Timept, y = m1, color = Genotype)) + geom_point() + stat_summary(fun.data=mean_se, fun = mean, geom = 'line', aes(group = Genotype)) + theme_minimal(base_size = 20)
# Overlay m1 model over summary data
ggplot(d3, aes(Timept, log10bact, color=Genotype)) +
stat_summary(fun.data=mean_se, geom="pointrange") +
stat_summary(aes(y=m1, linetype=Genotype), fun=mean, geom="line")+theme_minimal(base_size = 20)
# show info about m2 model fit
summary(poly_m.2)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: log10bact ~ (ot1 + ot2) + Genotype + ot1:Genotype + ot2:Genotype +
## (1 + ot1 + ot2 | TechRep)
## Data: d3
##
## AIC BIC logLik deviance df.resid
## 114.3 142.4 -44.1 88.3 51
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.4035 -0.5224 0.1970 0.6464 1.3560
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## TechRep (Intercept) 0.000e+00 0.000e+00
## ot1 1.508e-11 3.883e-06 NaN
## ot2 1.116e-10 1.056e-05 NaN -0.93
## Residual 2.326e-01 4.823e-01
## Number of obs: 64, groups: TechRep, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 5.15229 0.08526 60.429
## ot1 1.00421 0.24116 4.164
## ot2 0.56772 0.24116 2.354
## GenotypeWT -0.06476 0.12058 -0.537
## ot1:GenotypeWT -0.04431 0.34105 -0.130
## ot2:GenotypeWT -0.08814 0.34105 -0.258
##
## Correlation of Fixed Effects:
## (Intr) ot1 ot2 GntyWT o1:GWT
## ot1 0.000
## ot2 0.000 0.000
## GenotypeWT -0.707 0.000 0.000
## ot1:GntypWT 0.000 -0.707 0.000 0.000
## ot2:GntypWT 0.000 0.000 -0.707 0.000 0.000
## convergence code: 0
## boundary (singular) fit: see ?isSingular
# get parameter estimates and estimate p-values
m.coefs <- data.frame(coef(summary(poly_m.2)))
m.coefs$p <- 2*(1-pnorm(abs(m.coefs$t.value)))
m.coefs
## Estimate Std..Error t.value p
## (Intercept) 5.15229031 0.0852614 60.4293390 0.000000e+00
## ot1 1.00421198 0.2411557 4.1641649 3.124942e-05
## ot2 0.56772021 0.2411557 2.3541649 1.856437e-02
## GenotypeWT -0.06476083 0.1205778 -0.5370873 5.912073e-01
## ot1:GenotypeWT -0.04431146 0.3410456 -0.1299283 8.966232e-01
## ot2:GenotypeWT -0.08814389 0.3410456 -0.2584519 7.960582e-01
# Plot m2 model predictions
ggplot(d3, aes(x = Timept, y = m2, color = Genotype)) + geom_point() + stat_summary(fun.data=mean_se, fun = mean, geom = 'line', aes(group = Genotype)) + theme_minimal(base_size = 20)
# Overlay m2 model over summary data
ggplot(d3, aes(Timept, log10bact, color=Genotype)) +
stat_summary(fun.data=mean_se, geom="pointrange") +
stat_summary(aes(y=m2, linetype=Genotype), fun=mean, geom="line")+theme_minimal(base_size = 20)
##### GCA - Cubic ######
d3 <- d2_facet
# add timeBin
orderedTimepts <- unique(d3$Timept)[order(unique(d3$Timept))]
d3$timeBin <- match(d3$Timept,orderedTimepts)
# add 3rd degree orthogonal polynomial
t <- poly(unique(d3$Timept), 3) # one row for each timepoint, three cols for three dimensions
d3[,paste("ot", 1:3, sep="")] <- t[d3$timeBin, 1:3] # placeholders in ot1, ot2, ot3 for now
# Create models
# no group effects model
# poly_m.base <- lmer(log10bact ~ (ot1+ot2+ot3) + (1+ot1+ot2+ot3|BiolRepl/TechRep), data=d3, REML = FALSE)
poly_m.base <- lmer(log10bact ~ (ot1+ot2+ot3) + (1+ot1+ot2+ot3|TechRep), data=d3, REML = FALSE)
## boundary (singular) fit: see ?isSingular
# add effect of Genotype on intercept
# poly_m.0 <- lmer(log10bact ~ (ot1+ot2+ot3) + Genotype + (1+ot1+ot2+ot3|BiolRepl/TechRep), data=d3, REML = FALSE)
poly_m.0 <- lmer(log10bact ~ (ot1+ot2+ot3) + Genotype + (1+ot1+ot2+ot3|TechRep), data=d3, REML = FALSE)
## boundary (singular) fit: see ?isSingular
# add effect of Genotype on linear term
# poly_m.1 <- lmer(log10bact ~ (ot1+ot2+ot3) + Genotype + ot1:Genotype + (1+ot1+ot2+ot3|BiolRepl/TechRep), data=d3, REML = FALSE)
poly_m.1 <- lmer(log10bact ~ (ot1+ot2+ot3) + Genotype + ot1:Genotype + (1+ot1+ot2+ot3|TechRep), data=d3, REML = FALSE)
## boundary (singular) fit: see ?isSingular
# add effect of Genotype on quadratic term
# poly_m.2 <- lmer(log10bact ~ (ot1+ot2+ot3) + Genotype + ot1:Genotype +ot2:Genotype + (1+ot1+ot2+ot3|BiolRepl/TechRep), data=d3, REML = FALSE)
poly_m.2 <- lmer(log10bact ~ (ot1+ot2+ot3) + Genotype + ot1:Genotype +ot2:Genotype + (1+ot1+ot2+ot3|TechRep), data=d3, REML = FALSE)
## boundary (singular) fit: see ?isSingular
# add effect of Genotype on cubic term
# poly_m.3 <- lmer(log10bact ~ (ot1+ot2+ot3) + (ot1+ot2+ot3)*Genotype + (1+ot1+ot2+ot3|BiolRepl/TechRep), data=d3, REML = FALSE)
poly_m.3 <- lmer(log10bact ~ (ot1+ot2+ot3) + (ot1+ot2+ot3)*Genotype + (1+ot1+ot2+ot3|TechRep), data=d3, REML = FALSE)
## boundary (singular) fit: see ?isSingular
# compare the goodness of fit of the models
anova(poly_m.base, poly_m.0, poly_m.1, poly_m.2, poly_m.3)
## Data: d3
## Models:
## poly_m.base: log10bact ~ (ot1 + ot2 + ot3) + (1 + ot1 + ot2 + ot3 | TechRep)
## poly_m.0: log10bact ~ (ot1 + ot2 + ot3) + Genotype + (1 + ot1 + ot2 + ot3 |
## poly_m.0: TechRep)
## poly_m.1: log10bact ~ (ot1 + ot2 + ot3) + Genotype + ot1:Genotype + (1 +
## poly_m.1: ot1 + ot2 + ot3 | TechRep)
## poly_m.2: log10bact ~ (ot1 + ot2 + ot3) + Genotype + ot1:Genotype + ot2:Genotype +
## poly_m.2: (1 + ot1 + ot2 + ot3 | TechRep)
## poly_m.3: log10bact ~ (ot1 + ot2 + ot3) + (ot1 + ot2 + ot3) * Genotype +
## poly_m.3: (1 + ot1 + ot2 + ot3 | TechRep)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## poly_m.base 15 52.846 85.230 -11.423 22.846
## poly_m.0 16 54.039 88.581 -11.020 22.039 0.8071 1 0.3690
## poly_m.1 17 55.992 92.693 -10.996 21.992 0.0475 1 0.8274
## poly_m.2 18 57.803 96.663 -10.902 21.803 0.1885 1 0.6642
## poly_m.3 19 59.763 100.781 -10.881 21.763 0.0407 1 0.8402
#Add model fits to d3 dataframe
d3$base <- predict(poly_m.base)
d3$m0 <- predict(poly_m.0)
d3$m1 <- predict(poly_m.1)
d3$m2 <- predict(poly_m.2)
d3$m3 <- predict(poly_m.3)
## Show models on plots
# show info about base model fit
summary(poly_m.base)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: log10bact ~ (ot1 + ot2 + ot3) + (1 + ot1 + ot2 + ot3 | TechRep)
## Data: d3
##
## AIC BIC logLik deviance df.resid
## 52.8 85.2 -11.4 22.8 49
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.5404 -0.6480 0.3108 0.7951 1.2008
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## TechRep (Intercept) 0.000e+00 0.000e+00
## ot1 1.304e-12 1.142e-06 NaN
## ot2 4.948e-13 7.034e-07 NaN 0.68
## ot3 3.883e-13 6.232e-07 NaN 0.31 -0.42
## Residual 8.367e-02 2.893e-01
## Number of obs: 64, groups: TechRep, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 5.11991 0.03616 141.603
## ot1 0.98206 0.10227 9.603
## ot2 0.52365 0.10227 5.120
## ot3 -1.09657 0.10227 -10.723
##
## Correlation of Fixed Effects:
## (Intr) ot1 ot2
## ot1 0.000
## ot2 0.000 0.000
## ot3 0.000 0.000 0.000
## convergence code: 0
## boundary (singular) fit: see ?isSingular
# get parameter estimates and estimate p-values
m.coefs <- data.frame(coef(summary(poly_m.base)))
m.coefs$p <- 2*(1-pnorm(abs(m.coefs$t.value)))
m.coefs
## Estimate Std..Error t.value p
## (Intercept) 5.1199099 0.03615672 141.603261 0.000000e+00
## ot1 0.9820562 0.10226666 9.602898 0.000000e+00
## ot2 0.5236483 0.10226666 5.120420 3.048554e-07
## ot3 -1.0965733 0.10226666 -10.722686 0.000000e+00
# Plot base model predictions
ggplot(d3, aes(x = Timept, y = base, color = Genotype)) + geom_point() + stat_summary(fun.data=mean_se, fun = mean, geom = 'line', aes(group = Genotype)) + theme_minimal(base_size = 20)
# Overlay base model over summary data
ggplot(d3, aes(Timept, log10bact, color=Genotype)) +
stat_summary(fun.data=mean_se, geom="pointrange") +
stat_summary(aes(y=base, linetype=Genotype), fun=mean, geom="line")+theme_minimal(base_size = 20)
# show info about m0 model fit
summary(poly_m.0)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: log10bact ~ (ot1 + ot2 + ot3) + Genotype + (1 + ot1 + ot2 + ot3 |
## TechRep)
## Data: d3
##
## AIC BIC logLik deviance df.resid
## 54.0 88.6 -11.0 22.0 48
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.5568 -0.6667 0.3442 0.6983 1.2761
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## TechRep (Intercept) 0.000e+00 0.000e+00
## ot1 3.447e-12 1.857e-06 NaN
## ot2 7.836e-12 2.799e-06 NaN -0.70
## ot3 6.156e-12 2.481e-06 NaN -0.02 0.04
## Residual 8.262e-02 2.874e-01
## Number of obs: 64, groups: TechRep, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 5.15229 0.05081 101.399
## ot1 0.98206 0.10162 9.664
## ot2 0.52365 0.10162 5.153
## ot3 -1.09657 0.10162 -10.791
## GenotypeWT -0.06476 0.07186 -0.901
##
## Correlation of Fixed Effects:
## (Intr) ot1 ot2 ot3
## ot1 0.000
## ot2 0.000 0.000
## ot3 0.000 0.000 0.000
## GenotypeWT -0.707 0.000 0.000 0.000
## convergence code: 0
## boundary (singular) fit: see ?isSingular
# get parameter estimates and estimate p-values
m.coefs <- data.frame(coef(summary(poly_m.0)))
m.coefs$p <- 2*(1-pnorm(abs(m.coefs$t.value)))
m.coefs
## Estimate Std..Error t.value p
## (Intercept) 5.15229031 0.05081193 101.3992311 0.000000e+00
## ot1 0.98205625 0.10162385 9.6636391 0.000000e+00
## ot2 0.52364827 0.10162385 5.1528086 2.566138e-07
## ot3 -1.09657327 0.10162385 -10.7905106 0.000000e+00
## GenotypeWT -0.06476083 0.07185892 -0.9012219 3.674704e-01
# Plot m0 model predictions
ggplot(d3, aes(x = Timept, y = m0, color = Genotype)) + geom_point() + stat_summary(fun.data=mean_se, fun = mean, geom = 'line', aes(group = Genotype)) + theme_minimal(base_size = 20)
# Overlay m0 model over summary data
ggplot(d3, aes(Timept, log10bact, color=Genotype)) +
stat_summary(fun.data=mean_se, geom="pointrange") +
stat_summary(aes(y=m0, linetype=Genotype), fun=mean, geom="line")+theme_minimal(base_size = 20)
# show info about m1 model fit
summary(poly_m.1)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: log10bact ~ (ot1 + ot2 + ot3) + Genotype + ot1:Genotype + (1 +
## ot1 + ot2 + ot3 | TechRep)
## Data: d3
##
## AIC BIC logLik deviance df.resid
## 56.0 92.7 -11.0 22.0 47
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.5293 -0.6519 0.3679 0.7389 1.2713
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## TechRep (Intercept) 0.000e+00 0.000e+00
## ot1 5.968e-13 7.725e-07 NaN
## ot2 1.413e-11 3.759e-06 NaN 0.91
## ot3 3.549e-11 5.957e-06 NaN 0.40 0.31
## Residual 8.256e-02 2.873e-01
## Number of obs: 64, groups: TechRep, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 5.15229 0.05079 101.437
## ot1 1.00421 0.14366 6.990
## ot2 0.52365 0.10159 5.155
## ot3 -1.09657 0.10159 -10.795
## GenotypeWT -0.06476 0.07183 -0.902
## ot1:GenotypeWT -0.04431 0.20317 -0.218
##
## Correlation of Fixed Effects:
## (Intr) ot1 ot2 ot3 GntyWT
## ot1 0.000
## ot2 0.000 0.000
## ot3 0.000 0.000 0.000
## GenotypeWT -0.707 0.000 0.000 0.000
## ot1:GntypWT 0.000 -0.707 0.000 0.000 0.000
## convergence code: 0
## boundary (singular) fit: see ?isSingular
# get parameter estimates and estimate p-values
m.coefs <- data.frame(coef(summary(poly_m.1)))
m.coefs$p <- 2*(1-pnorm(abs(m.coefs$t.value)))
m.coefs
## Estimate Std..Error t.value p
## (Intercept) 5.15229031 0.05079305 101.4369056 0.000000e+00
## ot1 1.00421198 0.14366445 6.9899822 2.749134e-12
## ot2 0.52364827 0.10158611 5.1547231 2.540060e-07
## ot3 -1.09657327 0.10158611 -10.7945198 0.000000e+00
## GenotypeWT -0.06476083 0.07183223 -0.9015567 3.672924e-01
## ot1:GenotypeWT -0.04431146 0.20317222 -0.2180980 8.273527e-01
# Plot m1 model predictions
ggplot(d3, aes(x = Timept, y = m1, color = Genotype)) + geom_point() + stat_summary(fun.data=mean_se, fun = mean, geom = 'line', aes(group = Genotype)) + theme_minimal(base_size = 20)
# Overlay m1 model over summary data
ggplot(d3, aes(Timept, log10bact, color=Genotype)) +
stat_summary(fun.data=mean_se, geom="pointrange") +
stat_summary(aes(y=m1, linetype=Genotype), fun=mean, geom="line")+theme_minimal(base_size = 20)
# show info about m2 model fit
summary(poly_m.2)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula:
## log10bact ~ (ot1 + ot2 + ot3) + Genotype + ot1:Genotype + ot2:Genotype +
## (1 + ot1 + ot2 + ot3 | TechRep)
## Data: d3
##
## AIC BIC logLik deviance df.resid
## 57.8 96.7 -10.9 21.8 46
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.5336 -0.6726 0.3752 0.7076 1.2149
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## TechRep (Intercept) 0.000e+00 0.000e+00
## ot1 5.463e-14 2.337e-07 NaN
## ot2 1.789e-11 4.230e-06 NaN 0.28
## ot3 2.067e-11 4.547e-06 NaN 0.93 0.17
## Residual 8.232e-02 2.869e-01
## Number of obs: 64, groups: TechRep, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 5.15229 0.05072 101.586
## ot1 1.00421 0.14345 7.000
## ot2 0.56772 0.14345 3.958
## ot3 -1.09657 0.10144 -10.810
## GenotypeWT -0.06476 0.07173 -0.903
## ot1:GenotypeWT -0.04431 0.20287 -0.218
## ot2:GenotypeWT -0.08814 0.20287 -0.434
##
## Correlation of Fixed Effects:
## (Intr) ot1 ot2 ot3 GntyWT o1:GWT
## ot1 0.000
## ot2 0.000 0.000
## ot3 0.000 0.000 0.000
## GenotypeWT -0.707 0.000 0.000 0.000
## ot1:GntypWT 0.000 -0.707 0.000 0.000 0.000
## ot2:GntypWT 0.000 0.000 -0.707 0.000 0.000 0.000
## convergence code: 0
## boundary (singular) fit: see ?isSingular
# get parameter estimates and estimate p-values
m.coefs <- data.frame(coef(summary(poly_m.2)))
m.coefs$p <- 2*(1-pnorm(abs(m.coefs$t.value)))
m.coefs
## Estimate Std..Error t.value p
## (Intercept) 5.15229031 0.05071831 101.5863918 0.000000e+00
## ot1 1.00421198 0.14345305 7.0002833 2.554401e-12
## ot2 0.56772021 0.14345305 3.9575333 7.572775e-05
## ot3 -1.09657327 0.10143662 -10.8104275 0.000000e+00
## GenotypeWT -0.06476083 0.07172652 -0.9028853 3.665868e-01
## ot1:GenotypeWT -0.04431146 0.20287325 -0.2184194 8.271023e-01
## ot2:GenotypeWT -0.08814389 0.20287325 -0.4344777 6.639416e-01
# Plot m2 model predictions
ggplot(d3, aes(x = Timept, y = m2, color = Genotype)) + geom_point() + stat_summary(fun.data=mean_se, fun = mean, geom = 'line', aes(group = Genotype)) + theme_minimal(base_size = 20)
# Overlay m2 model over summary data
ggplot(d3, aes(Timept, log10bact, color=Genotype)) +
stat_summary(fun.data=mean_se, geom="pointrange") +
stat_summary(aes(y=m2, linetype=Genotype), fun=mean, geom="line")+theme_minimal(base_size = 20)
# show info about m3 model fit
summary(poly_m.3)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: log10bact ~ (ot1 + ot2 + ot3) + (ot1 + ot2 + ot3) * Genotype +
## (1 + ot1 + ot2 + ot3 | TechRep)
## Data: d3
##
## AIC BIC logLik deviance df.resid
## 59.8 100.8 -10.9 21.8 45
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.5063 -0.6448 0.3829 0.7368 1.2069
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## TechRep (Intercept) 0.000e+00 0.000e+00
## ot1 1.071e-11 3.272e-06 NaN
## ot2 8.643e-12 2.940e-06 NaN 0.99
## ot3 5.329e-11 7.300e-06 NaN 0.51 0.49
## Residual 8.226e-02 2.868e-01
## Number of obs: 64, groups: TechRep, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 5.15229 0.05070 101.619
## ot1 1.00421 0.14341 7.003
## ot2 0.56772 0.14341 3.959
## ot3 -1.11702 0.14341 -7.789
## GenotypeWT -0.06476 0.07170 -0.903
## ot1:GenotypeWT -0.04431 0.20281 -0.218
## ot2:GenotypeWT -0.08814 0.20281 -0.435
## ot3:GenotypeWT 0.04090 0.20281 0.202
##
## Correlation of Fixed Effects:
## (Intr) ot1 ot2 ot3 GntyWT o1:GWT o2:GWT
## ot1 0.000
## ot2 0.000 0.000
## ot3 0.000 0.000 0.000
## GenotypeWT -0.707 0.000 0.000 0.000
## ot1:GntypWT 0.000 -0.707 0.000 0.000 0.000
## ot2:GntypWT 0.000 0.000 -0.707 0.000 0.000 0.000
## ot3:GntypWT 0.000 0.000 0.000 -0.707 0.000 0.000 0.000
## convergence code: 0
## boundary (singular) fit: see ?isSingular
# get parameter estimates and estimate p-values
m.coefs <- data.frame(coef(summary(poly_m.3)))
m.coefs$p <- 2*(1-pnorm(abs(m.coefs$t.value)))
m.coefs
## Estimate Std..Error t.value p
## (Intercept) 5.15229031 0.05070220 101.6186661 0.000000e+00
## ot1 1.00421198 0.14340749 7.0025073 2.514211e-12
## ot2 0.56772021 0.14340749 3.9587906 7.533025e-05
## ot3 -1.11702390 0.14340749 -7.7891602 6.661338e-15
## GenotypeWT -0.06476083 0.07170374 -0.9031722 3.664345e-01
## ot1:GenotypeWT -0.04431146 0.20280881 -0.2184888 8.270483e-01
## ot2:GenotypeWT -0.08814389 0.20280881 -0.4346157 6.638414e-01
## ot3:GenotypeWT 0.04090126 0.20280881 0.2016740 8.401716e-01
# Plot m3 model predictions
ggplot(d3, aes(x = Timept, y = m3, color = Genotype)) + geom_point() + stat_summary(fun.data=mean_se, fun = mean, geom = 'line', aes(group = Genotype)) + theme_minimal(base_size = 20)
# Overlay m3 model over summary data
ggplot(d3, aes(Timept, log10bact, color=Genotype)) +
stat_summary(fun.data=mean_se, geom="pointrange") +
stat_summary(aes(y=m3, linetype=Genotype), fun=mean, geom="line")+theme_minimal(base_size = 20)
