Parameters
library(lme4)
## Loading required package: Matrix
library(lattice)
cnucl= c(0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.67, 0.7, 0.8, 0.9, 1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2)
seed=1:100
directoryname <- "20181201"
## Path to the directory
wd <- paste0("/homedir/lgay/sead/nico/", directoryname, "/outfiles")
setwd(wd)
Compile AICc for each seed/muSA combination
outfilename <- paste0("../outputdata", Sys.Date(), ".txt")
length(grep(".txt",list.files("./")))==length(cnucl)*length(seed)
outputdata=NULL
#seed_i=1
#cnucl_i=0
for (seed_i in seed){
for (cnucl_i in cnucl){
filename <- paste0("./simoutput_", directoryname, "seed", seed_i, "cnucl", cnucl_i,".txt")
if(file.exists(filename)){
outputdata <- read.table(file=filename, sep=" ", header=TRUE)
## Write data
if(file.exists(outfilename)){
write.table(outputdata, file=outfilename, append = TRUE, quote=FALSE, col.names=FALSE, row.names=FALSE, sep=",")
}else{
write.table(t(names(outputdata)), file=outfilename, quote=FALSE, col.names=FALSE, row.names=FALSE, sep=",")
write.table(outputdata, file=outfilename, append = TRUE, quote=FALSE, col.names=FALSE, row.names=FALSE, sep=",")
}
}else{
print(paste0("Output data for simulation with ", filename, " does not exist."))
}
}
}
Load data
#outfilename <- paste0("../outputdata", "2018-09-28", ".txt")
outfilename <- paste0("./outputdata", Sys.Date(), ".txt")
sim <- read.table(file=outfilename, sep=",", header=TRUE)
head(sim)
## seed cnuclsim ll_0.1_0 ll_0.2_0 ll_0.25_0 ll_0.3_0 ll_0.35_0 ll_0.4_0
## 1 1 0.1 2812.061 2812.061 2812.061 2812.061 2812.061 2812.061
## 2 1 0.2 2810.320 2810.320 2810.320 2810.320 2810.320 2810.320
## 3 1 0.3 2806.319 2806.074 2806.165 2806.139 2806.122 2806.113
## 4 1 0.4 2802.167 2802.677 2802.716 2802.646 2802.548 2802.455
## 5 1 0.5 2797.423 2798.577 2798.983 2799.235 2799.344 2799.342
## 6 1 0.6 2792.460 2794.006 2794.724 2795.324 2795.771 2796.061
## ll_0.45_0 ll_0.5_0 ll_0.55_0 ll_0.57_0 ll_0.59_0 ll_0.61_0 ll_0.63_0
## 1 2812.061 2812.061 2812.061 2812.061 2812.061 2812.061 2812.061
## 2 2810.320 2810.320 2810.320 2810.320 2810.320 2810.320 2810.320
## 3 2806.105 2806.106 2806.105 2806.105 2806.105 2806.105 2806.105
## 4 2802.372 2802.297 2802.232 2802.209 2802.186 2802.165 2802.145
## 5 2799.265 2799.160 2799.053 2799.012 2798.971 2798.931 2798.892
## 6 2796.214 2796.260 2796.230 2796.203 2796.169 2796.129 2796.088
## ll_0.65_0 ll_0.67_0 ll_0.69_0 ll_0.71_0 ll_0.73_0 ll_0.75_0 ll_0.77_0
## 1 2812.061 2812.061 2812.061 2812.061 2812.061 2812.061 2812.061
## 2 2810.320 2810.320 2810.320 2810.320 2810.320 2810.320 2810.320
## 3 2806.105 2806.105 2806.105 2806.105 2806.105 2806.105 2806.105
## 4 2802.126 2802.108 2802.091 2802.074 2802.059 2802.044 2802.030
## 5 2798.855 2798.819 2798.784 2798.750 2798.717 2798.685 2798.655
## 6 2796.045 2796.002 2795.959 2795.916 2795.873 2795.831 2795.790
## ll_0.79_0 ll_0.85_0 ll_0.9_0 ll_0.95_0 ll_1_0 ll_1.05_0 ll_1.1_0
## 1 2812.061 2812.061 2812.061 2812.061 2812.061 2812.061 2812.061
## 2 2810.320 2810.320 2810.320 2810.320 2810.320 2810.320 2810.320
## 3 2806.105 2806.105 2806.105 2806.105 2805.904 2806.105 2806.105
## 4 2802.016 2801.979 2801.952 2801.927 2801.905 2801.885 2801.867
## 5 2798.625 2798.543 2798.481 2798.425 2798.373 2798.325 2798.280
## 6 2795.749 2795.633 2795.542 2795.458 2795.379 2795.306 2795.237
## ll_1.15_0 ll_1.2_0 ll_1.25_0 ll_1.3_0 ll_1.4_0 ll_1.5_0 X.Intercept.
## 1 2812.061 2812.061 2812.061 2812.061 2812.061 2812.061 1.918722
## 2 2810.320 2810.320 2810.320 2810.320 2810.320 2810.320 1.933986
## 3 2806.105 2805.904 2806.105 2806.105 2806.105 2806.105 1.955839
## 4 2801.851 2801.835 2801.821 2801.809 2801.786 2801.766 1.969211
## 5 2798.240 2798.202 2798.167 2798.134 2798.075 2798.023 1.986153
## 6 2795.173 2795.114 2795.058 2795.005 2794.910 2794.826 2.004913
## TypeHomokaryon DonorSG1 geneticdistancestd cnucl cmit
## 1 -0.2177824 -0.1975147 -0.08739660 1.50 0
## 2 -0.2193891 -0.1950676 -0.08280780 0.75 0
## 3 -0.2458016 -0.1910054 -0.08371218 0.10 0
## 4 -0.2191639 -0.1887379 -0.07162940 0.25 0
## 5 -0.2182878 -0.1867543 -0.06676887 0.35 0
## 6 -0.2248205 -0.1854437 -0.06580383 0.50 0
## geno_synthesis_HetStockplates_HetPreculture_assay_plate
## 1 0.01395689
## 2 0.01395367
## 3 0.01395049
## 4 0.01395032
## 5 0.01395019
## 6 0.01395053
## geno_synthesis_HetStockplates_HetPreculture_assay geno
## 1 0.04657814 5.394631e-15
## 2 0.04695664 1.972965e-16
## 3 0.04639611 0.000000e+00
## 4 0.04636601 0.000000e+00
## 5 0.04635994 0.000000e+00
## 6 0.04634957 0.000000e+00
## Donor Donor.1 Midparent Acceptor Acceptormito Acceptormito.1
## 1 7.100408e-02 0.000000000 0.00000000 0 0.3257154 0.0000000000
## 2 7.069973e-02 0.000000000 0.00000000 0 0.3257160 0.0000000000
## 3 2.322946e-08 0.002854711 0.08629677 0 0.2678054 0.0009739218
## 4 6.167674e-09 0.004262198 0.10066582 0 0.2241363 0.0013862757
## 5 5.631444e-10 0.005262540 0.11744053 0 0.2098663 0.0016648742
## 6 1.139982e-10 0.005662783 0.11637131 0 0.2187236 0.0013338196
## Residual
## 1 0.01107095
## 2 0.01107095
## 3 0.01107096
## 4 0.01107096
## 5 0.01107096
## 6 0.01107096
## Section cluster specific
##Dataset for CI
simCI <- sim[sim$cnuclsim==0.67,]
head(simCI)
## seed cnuclsim ll_0.1_0 ll_0.2_0 ll_0.25_0 ll_0.3_0 ll_0.35_0 ll_0.4_0
## 7 1 0.67 2789.102 2790.658 2791.526 2792.335 2793.023 2793.550
## 28 2 0.67 2871.071 2875.499 2877.715 2879.783 2881.202 2881.894
## 68 4 0.67 2909.102 2910.568 2911.249 2911.880 2912.460 2912.989
## 88 5 0.67 2771.021 2773.753 2775.061 2776.233 2777.135 2777.665
## 107 6 0.67 2837.995 2840.814 2841.838 2843.018 2844.050 2844.713
## 127 7 0.67 2850.553 2852.772 2853.656 2854.334 2854.797 2855.058
## ll_0.45_0 ll_0.5_0 ll_0.55_0 ll_0.57_0 ll_0.59_0 ll_0.61_0 ll_0.63_0
## 7 2793.911 2794.128 2794.231 2794.247 2794.252 2794.246 2794.232
## 28 2882.092 2881.973 2881.652 2881.485 2881.302 2881.105 2880.898
## 68 2913.463 2913.871 2914.201 2914.309 2914.404 2914.484 2914.550
## 88 2777.810 2777.646 2777.257 2777.055 2776.833 2776.593 2776.340
## 107 2845.110 2845.326 2845.423 2845.438 2845.443 2845.440 2845.429
## 127 2855.147 2855.100 2854.955 2854.878 2854.792 2854.699 2854.600
## ll_0.65_0 ll_0.67_0 ll_0.69_0 ll_0.71_0 ll_0.73_0 ll_0.75_0 ll_0.77_0
## 7 2794.211 2794.183 2794.150 2794.113 2794.072 2794.031 2793.989
## 28 2880.682 2880.459 2880.231 2880.007 2879.798 2879.601 2879.416
## 68 2914.602 2914.642 2914.670 2914.686 2914.692 2914.691 2914.685
## 88 2776.096 2775.871 2775.662 2775.469 2775.290 2775.122 2774.965
## 107 2845.413 2845.393 2845.362 2845.339 2845.308 2845.275 2845.240
## 127 2854.496 2854.391 2854.289 2854.180 2854.093 2853.999 2853.909
## ll_0.79_0 ll_0.85_0 ll_0.9_0 ll_0.95_0 ll_1_0 ll_1.05_0 ll_1.1_0
## 7 2793.947 2793.821 2793.719 2793.621 2793.529 2793.441 2793.359
## 28 2879.241 2878.769 2878.429 2878.127 2875.736 2877.615 2877.396
## 68 2914.675 2914.620 2914.554 2914.475 2914.389 2914.298 2914.205
## 88 2774.819 2774.429 2774.151 2773.908 2773.694 2773.503 2773.332
## 107 2845.204 2845.090 2844.991 2844.893 2844.796 2844.701 2844.609
## 127 2853.821 2853.577 2853.392 2853.223 2853.069 2852.927 2852.797
## ll_1.15_0 ll_1.2_0 ll_1.25_0 ll_1.3_0 ll_1.4_0 ll_1.5_0 X.Intercept.
## 7 2793.281 2793.208 2793.139 2793.075 2792.957 2792.851 2.019490
## 28 2877.197 2877.016 2876.850 2876.698 2876.427 2876.195 1.391460
## 68 2914.112 2914.020 2913.929 2913.841 2913.673 2913.518 2.153991
## 88 2773.179 2773.039 2772.914 2772.799 2772.595 2772.420 1.660702
## 107 2844.520 2844.435 2844.354 2844.276 2844.132 2844.001 2.298031
## 127 2852.678 2852.567 2852.464 2852.369 2852.198 2852.048 1.619981
## TypeHomokaryon DonorSG1 geneticdistancestd cnucl cmit
## 7 -0.2394890 -0.1849476 -0.06742100 0.59 0
## 28 -0.6164464 -0.1387300 -0.06631511 0.45 0
## 68 -0.7360434 -0.1072474 -0.07007662 0.73 0
## 88 0.4178379 -0.1866393 -0.03126606 0.45 0
## 107 -1.0238620 -0.1619739 -0.18737317 0.59 0
## 127 -1.0368517 -0.1609939 -0.10063944 0.45 0
## geno_synthesis_HetStockplates_HetPreculture_assay_plate
## 7 0.01395079
## 28 0.01158641
## 68 0.01176485
## 88 0.01197161
## 107 0.01256663
## 127 0.01175614
## geno_synthesis_HetStockplates_HetPreculture_assay geno
## 7 0.04633488 0.000000e+00
## 28 0.04754702 4.375610e-16
## 68 0.03733567 0.000000e+00
## 88 0.04510909 3.293986e-12
## 107 0.04485994 1.865659e-03
## 127 0.03936663 9.007747e-04
## Donor Donor.1 Midparent Acceptor Acceptormito
## 7 1.555819e-09 5.828402e-03 0.1187831 0.000000000 0.223212160
## 28 4.319239e-14 5.112036e-15 0.3986335 0.000000000 0.205622680
## 68 1.719569e-03 6.470246e-03 0.1396303 0.007487128 0.135847566
## 88 0.000000e+00 0.000000e+00 0.2438091 0.001919413 0.137371302
## 107 1.450992e-01 0.000000e+00 0.1512735 0.001385401 0.003140341
## 127 0.000000e+00 9.522539e-04 0.1948833 0.010453940 0.164196233
## Acceptormito.1 Residual
## 7 1.130465e-03 0.01107096
## 28 2.994397e-03 0.01093918
## 68 7.864477e-03 0.01098093
## 88 7.337478e-03 0.01149024
## 107 5.790726e-04 0.01101054
## 127 1.043270e-10 0.01124586
simCI <- data.frame(seed=paste0("sim", simCI$seed), simCI[, (grep("X.Intercept.", colnames(simCI))):ncol(simCI)])
## Dataset for power
simAIC <- sim[sim$cnuclsim!=0.67,]
simAIC <- data.frame(seed=paste0("sim", simAIC$seed), cnuclsim=simAIC$cnuclsim, cnucl=simAIC$cnucl, cmit=simAIC$cmit, simAIC[, grep("ll_", colnames(simAIC))])
head(simAIC)
## seed cnuclsim cnucl cmit ll_0.1_0 ll_0.2_0 ll_0.25_0 ll_0.3_0 ll_0.35_0
## 1 sim1 0.1 1.50 0 2812.061 2812.061 2812.061 2812.061 2812.061
## 2 sim1 0.2 0.75 0 2810.320 2810.320 2810.320 2810.320 2810.320
## 3 sim1 0.3 0.10 0 2806.319 2806.074 2806.165 2806.139 2806.122
## 4 sim1 0.4 0.25 0 2802.167 2802.677 2802.716 2802.646 2802.548
## 5 sim1 0.5 0.35 0 2797.423 2798.577 2798.983 2799.235 2799.344
## 6 sim1 0.6 0.50 0 2792.460 2794.006 2794.724 2795.324 2795.771
## ll_0.4_0 ll_0.45_0 ll_0.5_0 ll_0.55_0 ll_0.57_0 ll_0.59_0 ll_0.61_0
## 1 2812.061 2812.061 2812.061 2812.061 2812.061 2812.061 2812.061
## 2 2810.320 2810.320 2810.320 2810.320 2810.320 2810.320 2810.320
## 3 2806.113 2806.105 2806.106 2806.105 2806.105 2806.105 2806.105
## 4 2802.455 2802.372 2802.297 2802.232 2802.209 2802.186 2802.165
## 5 2799.342 2799.265 2799.160 2799.053 2799.012 2798.971 2798.931
## 6 2796.061 2796.214 2796.260 2796.230 2796.203 2796.169 2796.129
## ll_0.63_0 ll_0.65_0 ll_0.67_0 ll_0.69_0 ll_0.71_0 ll_0.73_0 ll_0.75_0
## 1 2812.061 2812.061 2812.061 2812.061 2812.061 2812.061 2812.061
## 2 2810.320 2810.320 2810.320 2810.320 2810.320 2810.320 2810.320
## 3 2806.105 2806.105 2806.105 2806.105 2806.105 2806.105 2806.105
## 4 2802.145 2802.126 2802.108 2802.091 2802.074 2802.059 2802.044
## 5 2798.892 2798.855 2798.819 2798.784 2798.750 2798.717 2798.685
## 6 2796.088 2796.045 2796.002 2795.959 2795.916 2795.873 2795.831
## ll_0.77_0 ll_0.79_0 ll_0.85_0 ll_0.9_0 ll_0.95_0 ll_1_0 ll_1.05_0
## 1 2812.061 2812.061 2812.061 2812.061 2812.061 2812.061 2812.061
## 2 2810.320 2810.320 2810.320 2810.320 2810.320 2810.320 2810.320
## 3 2806.105 2806.105 2806.105 2806.105 2806.105 2805.904 2806.105
## 4 2802.030 2802.016 2801.979 2801.952 2801.927 2801.905 2801.885
## 5 2798.655 2798.625 2798.543 2798.481 2798.425 2798.373 2798.325
## 6 2795.790 2795.749 2795.633 2795.542 2795.458 2795.379 2795.306
## ll_1.1_0 ll_1.15_0 ll_1.2_0 ll_1.25_0 ll_1.3_0 ll_1.4_0 ll_1.5_0
## 1 2812.061 2812.061 2812.061 2812.061 2812.061 2812.061 2812.061
## 2 2810.320 2810.320 2810.320 2810.320 2810.320 2810.320 2810.320
## 3 2806.105 2806.105 2805.904 2806.105 2806.105 2806.105 2806.105
## 4 2801.867 2801.851 2801.835 2801.821 2801.809 2801.786 2801.766
## 5 2798.280 2798.240 2798.202 2798.167 2798.134 2798.075 2798.023
## 6 2795.237 2795.173 2795.114 2795.058 2795.005 2794.910 2794.826
load("./mbest.RData", verbose = FALSE)
Power
vpars <- function(m) {
nrow(m)*(nrow(m)+1)/2
}
llbest <- logLik(mbest)
csnucl <- attr(mbest, "csnucl")
csmit <- attr(mbest, "csmit")
csnucldf <- ifelse(csnucl!=1&csnucl!=0,1,0)
csmitdf <- ifelse(csmit!=1&csmit!=0,1,0)
ParNum <- sum(sapply(VarCorr(mbest), vpars))+length(fixef(mbest))+1+csnucldf+csmitdf
ParNumcnucl1 <- ParNum-1
numobs <- nrow(model.frame(mbest))
## Extract log-likelihood
simAIC$llmax <- apply(simAIC[, grep("ll_", colnames(simAIC))], 1, max)
## Compute AICc for best model and model with cnucl=1
simAIC$AICcmax <- (-2*simAIC$llmax+2*ParNum+(2*ParNum*(ParNum+1))/(numobs-ParNum-1))
simAIC$AICccnucl1 <- (-2*simAIC$ll_1_0+2*ParNumcnucl1+(2*ParNumcnucl1*(ParNumcnucl1+1))/(numobs-ParNumcnucl1-1))
simAIC$minAICc <- apply(simAIC[, c("AICcmax", "AICccnucl1")], 1, min)
simAIC$DeltaAICc <- simAIC$AICccnucl1 - simAIC$minAICc
## Compute statistical power
simAIC$DeltaAICcsup2 <- ifelse(simAIC$DeltaAICc > 2, 1, 0)
statpower <- tapply(simAIC$DeltaAICcsup2, simAIC$cnuclsim, mean)
pdf("./FigureS8_StatisticalPowerMGR.pdf")
plot(statpower~seq(0.1, 2, by=0.1), ylab="Statistical Power", xlab="cnucl", type="l", lwd=2, las=1, ylim=c(0, 1))
abline(h=0.8, lty=2)
dev.off()
## png
## 2
histogram(~DeltaAICc|cnuclsim, data=simAIC)

95% confidence interval
head(simCI)
## seed X.Intercept. TypeHomokaryon DonorSG1 geneticdistancestd cnucl
## 7 sim1 2.019490 -0.2394890 -0.1849476 -0.06742100 0.59
## 28 sim2 1.391460 -0.6164464 -0.1387300 -0.06631511 0.45
## 68 sim4 2.153991 -0.7360434 -0.1072474 -0.07007662 0.73
## 88 sim5 1.660702 0.4178379 -0.1866393 -0.03126606 0.45
## 107 sim6 2.298031 -1.0238620 -0.1619739 -0.18737317 0.59
## 127 sim7 1.619981 -1.0368517 -0.1609939 -0.10063944 0.45
## cmit geno_synthesis_HetStockplates_HetPreculture_assay_plate
## 7 0 0.01395079
## 28 0 0.01158641
## 68 0 0.01176485
## 88 0 0.01197161
## 107 0 0.01256663
## 127 0 0.01175614
## geno_synthesis_HetStockplates_HetPreculture_assay geno
## 7 0.04633488 0.000000e+00
## 28 0.04754702 4.375610e-16
## 68 0.03733567 0.000000e+00
## 88 0.04510909 3.293986e-12
## 107 0.04485994 1.865659e-03
## 127 0.03936663 9.007747e-04
## Donor Donor.1 Midparent Acceptor Acceptormito
## 7 1.555819e-09 5.828402e-03 0.1187831 0.000000000 0.223212160
## 28 4.319239e-14 5.112036e-15 0.3986335 0.000000000 0.205622680
## 68 1.719569e-03 6.470246e-03 0.1396303 0.007487128 0.135847566
## 88 0.000000e+00 0.000000e+00 0.2438091 0.001919413 0.137371302
## 107 1.450992e-01 0.000000e+00 0.1512735 0.001385401 0.003140341
## 127 0.000000e+00 9.522539e-04 0.1948833 0.010453940 0.164196233
## Acceptormito.1 Residual
## 7 1.130465e-03 0.01107096
## 28 2.994397e-03 0.01093918
## 68 7.864477e-03 0.01098093
## 88 7.337478e-03 0.01149024
## 107 5.790726e-04 0.01101054
## 127 1.043270e-10 0.01124586
## Compute acceptor, donor, heterokaryon and homokaryon variances
simCI$heterokaryonVar <- simCI$Donor.1 + simCI$Acceptor + ((simCI$cnucl^2)/2)*simCI$Midparent + simCI$Acceptormito.1 + simCI$geno
simCI$homokaryonVar <- simCI$Donor + simCI$Midparent + simCI$Acceptormito
simCI$DonorVar <- simCI$Donor.1 + simCI$Midparent*(simCI$cnucl/2)^2
simCI$AcceptorVar <- simCI$Acceptor + simCI$Midparent*(simCI$cnucl/2)^2 + simCI$Acceptormito.1
## Compute correlation between acceptor and donor genetic effects
simCI$AccDonCor <- (simCI$Midparent*(simCI$cnucl/2)^2) / sqrt(simCI$DonorVar*simCI$AcceptorVar)
simCI$homhetCor <- (simCI$Midparent*(simCI$cnucl/2))/ sqrt(simCI$heterokaryonVar*simCI$homokaryonVar)
## Slope of the regression of heterokaryon over mid-homokaryon parent value
simCI$midhomparhetslopeall <- (simCI$cnucl*simCI$Midparent)/(simCI$Midparent+simCI$Donor+simCI$Acceptormito+simCI$geno_synthesis_HetStockplates_HetPreculture_assay+simCI$geno_synthesis_HetStockplates_HetPreculture_assay_plate+simCI$Residual)
simCI$midhomparhetslopeaverage <- (simCI$cnucl*simCI$Midparent)/(simCI$Midparent+simCI$Donor+simCI$Acceptormito)
## The slope with the acceptor and donor homokaryon parents are equal as the mitochondrial loci have homokaryon- or heterokaryon-specific effects.
simCI$homparhetslopeall <- ((simCI$cnucl/2)*simCI$Midparent)/(simCI$Midparent+simCI$Donor+simCI$Acceptormito+simCI$geno_synthesis_HetStockplates_HetPreculture_assay+simCI$geno_synthesis_HetStockplates_HetPreculture_assay_plate+simCI$Residual)
simCI$homparhetslopeaverage <- ((simCI$cnucl/2)*simCI$Midparent)/(simCI$Midparent+simCI$Donor+simCI$Acceptormito)
head(simCI)
## seed X.Intercept. TypeHomokaryon DonorSG1 geneticdistancestd cnucl
## 7 sim1 2.019490 -0.2394890 -0.1849476 -0.06742100 0.59
## 28 sim2 1.391460 -0.6164464 -0.1387300 -0.06631511 0.45
## 68 sim4 2.153991 -0.7360434 -0.1072474 -0.07007662 0.73
## 88 sim5 1.660702 0.4178379 -0.1866393 -0.03126606 0.45
## 107 sim6 2.298031 -1.0238620 -0.1619739 -0.18737317 0.59
## 127 sim7 1.619981 -1.0368517 -0.1609939 -0.10063944 0.45
## cmit geno_synthesis_HetStockplates_HetPreculture_assay_plate
## 7 0 0.01395079
## 28 0 0.01158641
## 68 0 0.01176485
## 88 0 0.01197161
## 107 0 0.01256663
## 127 0 0.01175614
## geno_synthesis_HetStockplates_HetPreculture_assay geno
## 7 0.04633488 0.000000e+00
## 28 0.04754702 4.375610e-16
## 68 0.03733567 0.000000e+00
## 88 0.04510909 3.293986e-12
## 107 0.04485994 1.865659e-03
## 127 0.03936663 9.007747e-04
## Donor Donor.1 Midparent Acceptor Acceptormito
## 7 1.555819e-09 5.828402e-03 0.1187831 0.000000000 0.223212160
## 28 4.319239e-14 5.112036e-15 0.3986335 0.000000000 0.205622680
## 68 1.719569e-03 6.470246e-03 0.1396303 0.007487128 0.135847566
## 88 0.000000e+00 0.000000e+00 0.2438091 0.001919413 0.137371302
## 107 1.450992e-01 0.000000e+00 0.1512735 0.001385401 0.003140341
## 127 0.000000e+00 9.522539e-04 0.1948833 0.010453940 0.164196233
## Acceptormito.1 Residual heterokaryonVar homokaryonVar DonorVar
## 7 1.130465e-03 0.01107096 0.02763306 0.3419952 0.01616550
## 28 2.994397e-03 0.01093918 0.04335604 0.6042562 0.02018082
## 68 7.864477e-03 0.01098093 0.05902634 0.2771974 0.02507249
## 88 7.337478e-03 0.01149024 0.03394256 0.3811804 0.01234284
## 107 5.790726e-04 0.01101054 0.03015929 0.2995130 0.01316458
## 127 1.043270e-10 0.01124586 0.03203890 0.3590795 0.01081822
## AcceptorVar AccDonCor homhetCor midhomparhetslopeall
## 7 0.01146756 0.7592215 0.3604558 0.1695457
## 28 0.02317522 0.9331630 0.5541419 0.2660202
## 68 0.03395385 0.6375615 0.3984327 0.3022131
## 88 0.02159973 0.7559331 0.4822751 0.2439439
## 107 0.01512905 0.9328195 0.4695334 0.2425638
## 127 0.02031991 0.6654277 0.4088114 0.2080860
## midhomparhetslopeaverage homparhetslopeall homparhetslopeaverage
## 7 0.2049210 0.08477283 0.1024605
## 28 0.2968692 0.13301010 0.1484346
## 68 0.3677166 0.15110656 0.1838583
## 88 0.2878272 0.12197195 0.1439136
## 107 0.2979883 0.12128188 0.1489941
## 127 0.2442286 0.10404302 0.1221143
simboot <- simCI[, (grep("X.Intercept.", colnames(simCI))):ncol(simCI)]
head(simboot)
## X.Intercept. TypeHomokaryon DonorSG1 geneticdistancestd cnucl cmit
## 7 2.019490 -0.2394890 -0.1849476 -0.06742100 0.59 0
## 28 1.391460 -0.6164464 -0.1387300 -0.06631511 0.45 0
## 68 2.153991 -0.7360434 -0.1072474 -0.07007662 0.73 0
## 88 1.660702 0.4178379 -0.1866393 -0.03126606 0.45 0
## 107 2.298031 -1.0238620 -0.1619739 -0.18737317 0.59 0
## 127 1.619981 -1.0368517 -0.1609939 -0.10063944 0.45 0
## geno_synthesis_HetStockplates_HetPreculture_assay_plate
## 7 0.01395079
## 28 0.01158641
## 68 0.01176485
## 88 0.01197161
## 107 0.01256663
## 127 0.01175614
## geno_synthesis_HetStockplates_HetPreculture_assay geno
## 7 0.04633488 0.000000e+00
## 28 0.04754702 4.375610e-16
## 68 0.03733567 0.000000e+00
## 88 0.04510909 3.293986e-12
## 107 0.04485994 1.865659e-03
## 127 0.03936663 9.007747e-04
## Donor Donor.1 Midparent Acceptor Acceptormito
## 7 1.555819e-09 5.828402e-03 0.1187831 0.000000000 0.223212160
## 28 4.319239e-14 5.112036e-15 0.3986335 0.000000000 0.205622680
## 68 1.719569e-03 6.470246e-03 0.1396303 0.007487128 0.135847566
## 88 0.000000e+00 0.000000e+00 0.2438091 0.001919413 0.137371302
## 107 1.450992e-01 0.000000e+00 0.1512735 0.001385401 0.003140341
## 127 0.000000e+00 9.522539e-04 0.1948833 0.010453940 0.164196233
## Acceptormito.1 Residual heterokaryonVar homokaryonVar DonorVar
## 7 1.130465e-03 0.01107096 0.02763306 0.3419952 0.01616550
## 28 2.994397e-03 0.01093918 0.04335604 0.6042562 0.02018082
## 68 7.864477e-03 0.01098093 0.05902634 0.2771974 0.02507249
## 88 7.337478e-03 0.01149024 0.03394256 0.3811804 0.01234284
## 107 5.790726e-04 0.01101054 0.03015929 0.2995130 0.01316458
## 127 1.043270e-10 0.01124586 0.03203890 0.3590795 0.01081822
## AcceptorVar AccDonCor homhetCor midhomparhetslopeall
## 7 0.01146756 0.7592215 0.3604558 0.1695457
## 28 0.02317522 0.9331630 0.5541419 0.2660202
## 68 0.03395385 0.6375615 0.3984327 0.3022131
## 88 0.02159973 0.7559331 0.4822751 0.2439439
## 107 0.01512905 0.9328195 0.4695334 0.2425638
## 127 0.02031991 0.6654277 0.4088114 0.2080860
## midhomparhetslopeaverage homparhetslopeall homparhetslopeaverage
## 7 0.2049210 0.08477283 0.1024605
## 28 0.2968692 0.13301010 0.1484346
## 68 0.3677166 0.15110656 0.1838583
## 88 0.2878272 0.12197195 0.1439136
## 107 0.2979883 0.12128188 0.1489941
## 127 0.2442286 0.10404302 0.1221143
simboot <- simboot[,-c(which(colnames(simboot)=="heterokaryonVar"))]
bootmean <- apply(simboot, 2, mean)
bootsd <- data.frame(t(apply(simboot, 2, quantile, probs=c(0.025, 0.975), na.rm=TRUE)))
## Keep the slope separately from the rest
meanslope <- bootmean[(length(bootmean)-3):length(bootmean)]
slopeCI <- bootsd[(nrow(bootsd)-3):nrow(bootsd),]
Estimate <- bootmean[-((length(bootmean)-3):length(bootmean))]
bootsd <- bootsd[-((nrow(bootsd)-3):nrow(bootsd)),]
## Combine mean and CI
CI <- data.frame(Estimate=Estimate, bootsd)
CI <- CI[-c(which(rownames(CI)=="Donor"), which(rownames(CI)=="Donor.1"), which(rownames(CI)=="Acceptor"), which(rownames(CI)=="Midparent"), which(rownames(CI)=="cmit")),]
rownames(CI)[c(1:2, 4:3,5, 12, 14:13, 15:16, 8, 9:10, 7, 6, 11)]
## [1] "X.Intercept."
## [2] "TypeHomokaryon"
## [3] "geneticdistancestd"
## [4] "DonorSG1"
## [5] "cnucl"
## [6] "homokaryonVar"
## [7] "AcceptorVar"
## [8] "DonorVar"
## [9] "AccDonCor"
## [10] "homhetCor"
## [11] "geno"
## [12] "Acceptormito"
## [13] "Acceptormito.1"
## [14] "geno_synthesis_HetStockplates_HetPreculture_assay"
## [15] "geno_synthesis_HetStockplates_HetPreculture_assay_plate"
## [16] "Residual"
CI <- CI[c(1:2, 4:3,5, 12, 14:13, 15:16, 8, 9:10, 7, 6, 11),]
nam <- rep("", nrow(CI))
nam[1] <- "Fixed effects"
nam[nrow(CI)-10] <- "Random nuclear genetic effects"
nam[nrow(CI)-4] <- "Random mitochondrial genetic effects"
nam[nrow(CI)-2] <- "Random environmental effects"
Variable <- c("intercept","type (homokaryon)", "standardized genetic distance","donor slow grower","cnucl", "homokaryon", "acceptor","donor", "acceptor-donor correlation", "acceptor-homokaryon/donor-homokaryon correlation","acceptor x donor","homokaryon mitochondria","heterokaryon mitochondria", "assay", "plate","residual variance")
## Add model estimates
CI <- data.frame(Effect = nam, Variable, CI)
colnames(CI)[(ncol(CI)-1):ncol(CI)] <- c("2.5%CI", "97.5%CI")
write.csv(CI, file="./95%CIbestmodelGRstep1.csv", row.names=FALSE, quote=F)
##################################
## Proportion of variance explained
##################################
simCI2 <- simCI
head(simCI2)
## seed X.Intercept. TypeHomokaryon DonorSG1 geneticdistancestd cnucl
## 7 sim1 2.019490 -0.2394890 -0.1849476 -0.06742100 0.59
## 28 sim2 1.391460 -0.6164464 -0.1387300 -0.06631511 0.45
## 68 sim4 2.153991 -0.7360434 -0.1072474 -0.07007662 0.73
## 88 sim5 1.660702 0.4178379 -0.1866393 -0.03126606 0.45
## 107 sim6 2.298031 -1.0238620 -0.1619739 -0.18737317 0.59
## 127 sim7 1.619981 -1.0368517 -0.1609939 -0.10063944 0.45
## cmit geno_synthesis_HetStockplates_HetPreculture_assay_plate
## 7 0 0.01395079
## 28 0 0.01158641
## 68 0 0.01176485
## 88 0 0.01197161
## 107 0 0.01256663
## 127 0 0.01175614
## geno_synthesis_HetStockplates_HetPreculture_assay geno
## 7 0.04633488 0.000000e+00
## 28 0.04754702 4.375610e-16
## 68 0.03733567 0.000000e+00
## 88 0.04510909 3.293986e-12
## 107 0.04485994 1.865659e-03
## 127 0.03936663 9.007747e-04
## Donor Donor.1 Midparent Acceptor Acceptormito
## 7 1.555819e-09 5.828402e-03 0.1187831 0.000000000 0.223212160
## 28 4.319239e-14 5.112036e-15 0.3986335 0.000000000 0.205622680
## 68 1.719569e-03 6.470246e-03 0.1396303 0.007487128 0.135847566
## 88 0.000000e+00 0.000000e+00 0.2438091 0.001919413 0.137371302
## 107 1.450992e-01 0.000000e+00 0.1512735 0.001385401 0.003140341
## 127 0.000000e+00 9.522539e-04 0.1948833 0.010453940 0.164196233
## Acceptormito.1 Residual heterokaryonVar homokaryonVar DonorVar
## 7 1.130465e-03 0.01107096 0.02763306 0.3419952 0.01616550
## 28 2.994397e-03 0.01093918 0.04335604 0.6042562 0.02018082
## 68 7.864477e-03 0.01098093 0.05902634 0.2771974 0.02507249
## 88 7.337478e-03 0.01149024 0.03394256 0.3811804 0.01234284
## 107 5.790726e-04 0.01101054 0.03015929 0.2995130 0.01316458
## 127 1.043270e-10 0.01124586 0.03203890 0.3590795 0.01081822
## AcceptorVar AccDonCor homhetCor midhomparhetslopeall
## 7 0.01146756 0.7592215 0.3604558 0.1695457
## 28 0.02317522 0.9331630 0.5541419 0.2660202
## 68 0.03395385 0.6375615 0.3984327 0.3022131
## 88 0.02159973 0.7559331 0.4822751 0.2439439
## 107 0.01512905 0.9328195 0.4695334 0.2425638
## 127 0.02031991 0.6654277 0.4088114 0.2080860
## midhomparhetslopeaverage homparhetslopeall homparhetslopeaverage
## 7 0.2049210 0.08477283 0.1024605
## 28 0.2968692 0.13301010 0.1484346
## 68 0.3677166 0.15110656 0.1838583
## 88 0.2878272 0.12197195 0.1439136
## 107 0.2979883 0.12128188 0.1489941
## 127 0.2442286 0.10404302 0.1221143
ncol(simCI2)
## [1] 27
simCI2$HetPhenoVar <- simCI2$DonorVar + simCI2$AcceptorVar + simCI2$geno + simCI2$geno_synthesis_HetStockplates_HetPreculture_assay + simCI2$geno_synthesis_HetStockplates_HetPreculture_assay_plate + simCI2$Residual
simCI2$HomPhenoVar <- simCI2$homokaryonVar + simCI2$geno_synthesis_HetStockplates_HetPreculture_assay + simCI2$geno_synthesis_HetStockplates_HetPreculture_assay_plate + simCI2$Residual
simCI2$PropDon <- simCI2$DonorVar/simCI2$HetPhenoVar
simCI2$PropAcc <- (simCI2$AcceptorVar-simCI2$Acceptormito.1)/simCI2$HetPhenoVar
simCI2$Propgeno <- simCI2$geno/simCI2$HetPhenoVar
simCI2$PropAccMito <- simCI2$Acceptormito.1/simCI2$HetPhenoVar
simCI2$PropAssay <- simCI2$geno_synthesis_HetStockplates_HetPreculture_assay/simCI2$HetPhenoVar
simCI2$PropPlate <- simCI2$geno_synthesis_HetStockplates_HetPreculture_assay_plate/simCI2$HetPhenoVar
simCI2$PropRes <- simCI2$Residual/simCI2$HetPhenoVar
simCI2$PropHom<- (simCI2$homokaryonVar-simCI2$Acceptormito)/simCI2$HomPhenoVar
simCI2$PropHomMito <- simCI2$Acceptormito/simCI2$HomPhenoVar
simCI2$PropHomAssay <- simCI2$geno_synthesis_HetStockplates_HetPreculture_assay/simCI2$HomPhenoVar
simCI2$PropHomPlate <- simCI2$geno_synthesis_HetStockplates_HetPreculture_assay_plate/simCI2$HomPhenoVar
simCI2$PropHomRes <- simCI2$Residual/simCI2$HomPhenoVar
simCI2 <- simCI2[, grep("PropDon", (colnames(simCI2))):ncol(simCI2)]
bootmean <- apply(simCI2, 2, mean)
bootsd <- t(apply(simCI2, 2, quantile, probs=c(0.025, 0.975), na.rm=TRUE))
CI <- data.frame(Estimate=bootmean, bootsd)
names(Estimate)
## [1] "X.Intercept."
## [2] "TypeHomokaryon"
## [3] "DonorSG1"
## [4] "geneticdistancestd"
## [5] "cnucl"
## [6] "cmit"
## [7] "geno_synthesis_HetStockplates_HetPreculture_assay_plate"
## [8] "geno_synthesis_HetStockplates_HetPreculture_assay"
## [9] "geno"
## [10] "Donor"
## [11] "Donor.1"
## [12] "Midparent"
## [13] "Acceptor"
## [14] "Acceptormito"
## [15] "Acceptormito.1"
## [16] "Residual"
## [17] "homokaryonVar"
## [18] "DonorVar"
## [19] "AcceptorVar"
## [20] "AccDonCor"
## [21] "homhetCor"
rownames(bootsd)
## [1] "PropDon" "PropAcc" "Propgeno" "PropAccMito"
## [5] "PropAssay" "PropPlate" "PropRes" "PropHom"
## [9] "PropHomMito" "PropHomAssay" "PropHomPlate" "PropHomRes"
Variable <- c("acceptor nuclear","donor nuclear","acceptor x donor nuclear","heterokaryon mitochondria", "assay","plate","residual variance","homokaryon nuclear","homokaryon mitochondria", "assay","plate","residual variance")
nam <- rep("", nrow(CI))
nam[1] <- "Heterokaryon"
nam[nrow(CI)-4] <- "Homokaryon"
## Add model estimates
CI <- data.frame(Type=nam, Variable, CI)
colnames(CI)[(ncol(CI)-1):ncol(CI)] <- c("2.5%CI", "97.5%CI")
write.csv(CI, file="./Table4_95%CIPropVarGRstep1.csv", row.names=FALSE, quote=F)