1 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)

2 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."))
    }
  }
}

3 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)

4 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)

5 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)