1 Load script

library(MASS)
library(knitr)
source(purl(paste0("../Functions/","MMmodel.Rmd")))
## 
## 
## processing file: ../Functions/MMmodel.Rmd
## 
  |                                                                       
  |                                                                 |   0%
  |                                                                       
  |....                                                             |   6%
  |                                                                       
  |.......                                                          |  11%
  |                                                                       
  |...........                                                      |  17%
  |                                                                       
  |..............                                                   |  22%
  |                                                                       
  |..................                                               |  28%
  |                                                                       
  |......................                                           |  33%
  |                                                                       
  |.........................                                        |  39%
  |                                                                       
  |.............................                                    |  44%
  |                                                                       
  |................................                                 |  50%
  |                                                                       
  |....................................                             |  56%
  |                                                                       
  |........................................                         |  61%
  |                                                                       
  |...........................................                      |  67%
  |                                                                       
  |...............................................                  |  72%
  |                                                                       
  |...................................................              |  78%
  |                                                                       
  |......................................................           |  83%
  |                                                                       
  |..........................................................       |  89%
  |                                                                       
  |.............................................................    |  94%
  |                                                                       
  |.................................................................| 100%
## output file: MMmodel.R
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## The following object is masked from 'package:reshape':
## 
##     expand

dev.off()
## null device 
##           1
seedsim<- 1
#Nlocussim <- 5
Nlocussim <- 500
ngenosim <- 30
sdfitsim <- 0.001
nsimulCI <- 30
#stepc <- 0.1
stepc <- 0.025
start <- 0.1
stop <- 1
noautozygous <- FALSE

2 Constant h/s/c: Simulate data with h=0.25, s=0.01 and c=0.5

2.1 Prepare data

## Simulate data with constant haploid selection coeffcient and constant h
sim <- simdata(ngeno=ngenosim, nrep=10, Nlocus=Nlocussim, meanshap=0.01, shapesel=NULL, h=0.25, shapeh=NULL, c=0.5, shapec=NULL, mu=0.0001, sdfit=sdfitsim, seed=seedsim)
## [1] "Simulate data with seed=1"
## [1] "Simulations with constant shaploid, h and c across selected loci"
data <- sim$data
distmat <- sim$distmat
Expected_Hm00 <- sim$expectedH
datahap <- sim$datahap
datadip <- sim$datadip
Expected_Var <- sim$expectedVar
(observedc <- sim$observedc)
## [1] 0.5
simmat <- 1-distmat
datam00 <- data

2.2 Fit model

if(noautozygous){
  data<-data[!(data$Type=="diploid"&data$Acceptor==data$Donor),]
}

#m00 <- fitmidparentmodelcs(form="logfitness ~ 1  + distance + (1 | Midparent) + (0+  dummy(Type,\"diploid\")|Acceptor) + (0+  dummy(Type,\"diploid\") |Donor)+ (0+  dummy(Type,\"diploid\")|Acceptor_Donor) + (0+  dummy(Type,\"haploid\") |Donor)", cnucl=seq(0.5, 0.5, by=0.05), type="haploid", relmatrix=NA, REML=TRUE, data = data)

m00 <- fitmidparentmodelcs(form="logfitness ~ 1  + distance + (1 | Midparent) + (0+  dummy(Type,\"diploid\")|Acceptor) + (0+  dummy(Type,\"diploid\") |Donor)+ (0+  dummy(Type,\"diploid\")|Acceptor_Donor) + (0+  dummy(Type,\"haploid\") |Donor)", cnucl=seq(start, stop, by=stepc), type="haploid", relmatrix=NA, REML=TRUE, data = data)
## [1] "Fitting model with nuclei cs = 0.1"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## Warning in optwrap(optimizer, devfun, getStart(start, rho$lower, rho$pp), :
## convergence code 1 from bobyqa: bobyqa -- maximum number of function
## evaluations exceeded
## 'log Lik.' 26638.67 (df=8)
## [1] "Fitting model with nuclei cs = 0.125"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26723.97 (df=8)
## [1] "Fitting model with nuclei cs = 0.15"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26730.06 (df=8)
## [1] "Fitting model with nuclei cs = 0.175"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26648.52 (df=8)
## [1] "Fitting model with nuclei cs = 0.2"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26739.26 (df=8)
## [1] "Fitting model with nuclei cs = 0.225"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26742.81 (df=8)
## [1] "Fitting model with nuclei cs = 0.25"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26745.97 (df=8)
## [1] "Fitting model with nuclei cs = 0.275"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26748.99 (df=8)
## [1] "Fitting model with nuclei cs = 0.3"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26752.08 (df=8)
## [1] "Fitting model with nuclei cs = 0.325"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26755.43 (df=8)
## [1] "Fitting model with nuclei cs = 0.35"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26627.59 (df=8)
## [1] "Fitting model with nuclei cs = 0.375"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26633.46 (df=8)
## [1] "Fitting model with nuclei cs = 0.4"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26769.29 (df=8)
## [1] "Fitting model with nuclei cs = 0.425"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26647.9 (df=8)
## [1] "Fitting model with nuclei cs = 0.45"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26787.11 (df=8)
## [1] "Fitting model with nuclei cs = 0.475"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26804.86 (df=8)
## [1] "Fitting model with nuclei cs = 0.5"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26832.22 (df=8)
## [1] "Fitting model with nuclei cs = 0.525"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26675.43 (df=8)
## [1] "Fitting model with nuclei cs = 0.55"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26661.13 (df=8)
## [1] "Fitting model with nuclei cs = 0.575"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26650.89 (df=8)
## [1] "Fitting model with nuclei cs = 0.6"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26642.93 (df=8)
## [1] "Fitting model with nuclei cs = 0.625"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26750.62 (df=8)
## [1] "Fitting model with nuclei cs = 0.65"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26633.78 (df=8)
## [1] "Fitting model with nuclei cs = 0.675"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26737 (df=8)
## [1] "Fitting model with nuclei cs = 0.7"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26627.74 (df=8)
## [1] "Fitting model with nuclei cs = 0.725"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26625.42 (df=8)
## [1] "Fitting model with nuclei cs = 0.75"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26623.41 (df=8)
## [1] "Fitting model with nuclei cs = 0.775"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26621.66 (df=8)
## [1] "Fitting model with nuclei cs = 0.8"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26620.12 (df=8)
## [1] "Fitting model with nuclei cs = 0.825"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26618.87 (df=8)
## [1] "Fitting model with nuclei cs = 0.85"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26617.54 (df=8)
## [1] "Fitting model with nuclei cs = 0.875"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26616.51 (df=8)
## [1] "Fitting model with nuclei cs = 0.9"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26615.5 (df=8)
## [1] "Fitting model with nuclei cs = 0.925"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26614.57 (df=8)
## [1] "Fitting model with nuclei cs = 0.95"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26613.75 (df=8)
## [1] "Fitting model with nuclei cs = 0.975"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26613 (df=8)
## [1] "Fitting model with nuclei cs = 1"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26612.29 (df=8)
## [1] "Best nuclei c=0.5"
## [1] "95% confidence interval (0.44 0.51)"
#m00 <- fitmidparentmodelcs(form="logfitness ~ 1  + distance + (1 | Midparent) + (0+  dummy(Type,\"diploid\")|Acceptor) + (0+  dummy(Type,\"diploid\") |Donor)+ (0+  dummy(Type,\"diploid\")|Acceptor_Donor) + (0+  dummy(Type,\"haploid\") |Donor)", cnucl=seq(0.1, 0.9, by=0.1), type="haploid", relmatrix=list(Acceptor=simmat, Donor=simmat, Midparent=simmat), REML=TRUE, data = data)
#m00 <- fitmidparentmodelcs(form="logfitness ~ 1  + distance + (1 | Midparent) + (0+  dummy(Type,\"diploid\")|Acceptor) + (0+  dummy(Type,\"diploid\") |Donor)+ (0+  dummy(Type,\"diploid\")|Acceptor_Donor) + (0+  dummy(Type,\"haploid\") |Donor)", cnucl=0.5, type="haploid", relmatrix=list(Acceptor=simmat, Donor=simmat, Midparent=simmat), REML=TRUE, data = data)

summary(m00)
## Linear mixed model fit by REML ['lmerMod']
## 
## REML criterion at convergence: -53664.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5280 -0.6551 -0.0109  0.6801  3.5636 
## 
## Random effects:
##  Groups         Name                   Variance  Std.Dev. 
##  Acceptor_Donor dummy(Type, "diploid") 2.793e-09 5.285e-05
##  Donor          dummy(Type, "haploid") 6.044e-08 2.459e-04
##  Donor.1        dummy(Type, "diploid") 2.990e-20 1.729e-10
##  Acceptor       dummy(Type, "diploid") 0.000e+00 0.000e+00
##  Midparent      (Intercept)            9.331e-03 9.660e-02
##  Residual                              1.052e-06 1.026e-03
## Number of obs: 4950, groups:  
## Acceptor_Donor, 465; Donor, 30; Acceptor, 30; Midparent, 30
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept) 0.0004799  0.0001335   3.594
## distance    0.6255009  0.0015689 398.690
## 
## Correlation of Fixed Effects:
##          (Intr)
## distance -0.799

2.3 95%CI

sim <- SimCICompute(simdata=data, nsim=nsimulCI, mbest=m00)
CI <- data.frame(apply(sim, 2, quantile, probs=c(0.025, 0.975)))

write.table(CI, file=paste0('Constant_', Sys.Date(),".txt"), sep=" ", row.names = FALSE, quote=NA)

2.4 Check expected and fitted variance among haploid parents

## Check expected and fitted H
Fitted_Var <- data.frame(VarCorr(m00))$vcov[1]

c(Expected_Var=Expected_Var, Fitted_Var=Fitted_Var)
## Expected_Var   Fitted_Var 
##  8.71212e-04  2.79293e-09

2.5 Check expected and fitted H

## Check expected and fitted H
Fitted_Hm00=as.numeric(fixef(m00)[2])

c(Expected_H=Expected_Hm00, Fitted_H=Fitted_Hm00)
## Expected_H   Fitted_H 
##  0.6250000  0.6255009

2.6 Check expected and fitted cnucl

## Check expected and fitted cnucl
Fitted_cnuclm00 <- attr(m00, "csnucl") 
c(Expected_cnucl=0.5, Fitted_cnucl=Fitted_cnuclm00)
## Expected_cnucl   Fitted_cnucl 
##            0.5            0.5
## Compute observed c as ratio of autozygous diploid fitness over haploid fitness 
mean(datadip$meandiplogfitness/datahap$meanhaplogfitness)
## [1] 0.4972534
datadiphapratiom00 <- datadip$meandiplogfitness/datahap$meanhaplogfitness

## Likelihood profile
#attr(m00, "plotcnucl")

3 Constant h/variable s/constant c: Simulate data with distribution of selection coefficients across loci (mean s=0.01), but constant h=0.25, and c=0.5

3.1 Prepare data

sim <- simdata(ngeno=ngenosim, nrep=10, Nlocus=Nlocussim, meanshap=0.01, shapesel=2, h=0.25, shapeh=NULL, c=0.5, shapec=NULL, mu=0.0001, sdfit=sdfitsim, seed=seedsim)
## [1] "Simulate data with seed=1"
## [1] "Simulations with distribution of shaploid, h and/or c across selected loci"
## [1] "Mean of selection coefficients = 0.01"
## [1] "Standard deviation of selection coefficients = 0"
## [1] "Number of loci with p>1=1"
data <- sim$data
distmat <- sim$distmat
Expected_Hms <- sim$expectedH
datahap <- sim$datahap
datadip <- sim$datadip
Expected_Var <- sim$expectedVar
(observedc <- sim$observedc)
## [1] 0.5
simmat <- 1-distmat
datams <- data

3.2 Fit model

if(noautozygous){
  data<-data[!(data$Type=="diploid"&data$Acceptor==data$Donor),]
}


ms <- fitmidparentmodelcs(form="logfitness ~ 1  + distance + (1 | Midparent) + (0+  dummy(Type,\"diploid\")|Acceptor) + (0+  dummy(Type,\"diploid\") |Donor)+ (0+  dummy(Type,\"diploid\")|Acceptor_Donor) + (0+  dummy(Type,\"haploid\") |Donor)", cnucl=seq(start, stop, by=stepc), type="haploid", relmatrix=NA, REML=TRUE, data = data)
## [1] "Fitting model with nuclei cs = 0.1"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26028.9 (df=8)
## [1] "Fitting model with nuclei cs = 0.125"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26066.33 (df=8)
## [1] "Fitting model with nuclei cs = 0.15"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26072.96 (df=8)
## [1] "Fitting model with nuclei cs = 0.175"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26078.56 (df=8)
## [1] "Fitting model with nuclei cs = 0.2"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26083.38 (df=8)
## [1] "Fitting model with nuclei cs = 0.225"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26087.7 (df=8)
## [1] "Fitting model with nuclei cs = 0.25"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26091.71 (df=8)
## [1] "Fitting model with nuclei cs = 0.275"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26095.59 (df=8)
## [1] "Fitting model with nuclei cs = 0.3"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26099.35 (df=8)
## [1] "Fitting model with nuclei cs = 0.325"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26102.77 (df=8)
## [1] "Fitting model with nuclei cs = 0.35"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26105.37 (df=8)
## [1] "Fitting model with nuclei cs = 0.375"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26106.41 (df=8)
## [1] "Fitting model with nuclei cs = 0.4"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26105.36 (df=8)
## [1] "Fitting model with nuclei cs = 0.425"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26102.33 (df=8)
## [1] "Fitting model with nuclei cs = 0.45"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26098.01 (df=8)
## [1] "Fitting model with nuclei cs = 0.475"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26014.71 (df=8)
## [1] "Fitting model with nuclei cs = 0.5"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26011.43 (df=8)
## [1] "Fitting model with nuclei cs = 0.525"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26009.02 (df=8)
## [1] "Fitting model with nuclei cs = 0.55"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26078.33 (df=8)
## [1] "Fitting model with nuclei cs = 0.575"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26004.92 (df=8)
## [1] "Fitting model with nuclei cs = 0.6"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26003.5 (df=8)
## [1] "Fitting model with nuclei cs = 0.625"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26001.71 (df=8)
## [1] "Fitting model with nuclei cs = 0.65"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26079.18 (df=8)
## [1] "Fitting model with nuclei cs = 0.675"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 25999.22 (df=8)
## [1] "Fitting model with nuclei cs = 0.7"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 25998.16 (df=8)
## [1] "Fitting model with nuclei cs = 0.725"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 25997.18 (df=8)
## [1] "Fitting model with nuclei cs = 0.75"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 25996.32 (df=8)
## [1] "Fitting model with nuclei cs = 0.775"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 25995.75 (df=8)
## [1] "Fitting model with nuclei cs = 0.8"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 25994.87 (df=8)
## [1] "Fitting model with nuclei cs = 0.825"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26078.32 (df=8)
## [1] "Fitting model with nuclei cs = 0.85"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26078.23 (df=8)
## [1] "Fitting model with nuclei cs = 0.875"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26078.15 (df=8)
## [1] "Fitting model with nuclei cs = 0.9"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26078.07 (df=8)
## [1] "Fitting model with nuclei cs = 0.925"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26077.99 (df=8)
## [1] "Fitting model with nuclei cs = 0.95"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26077.92 (df=8)
## [1] "Fitting model with nuclei cs = 0.975"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26077.86 (df=8)
## [1] "Fitting model with nuclei cs = 1"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26077.79 (df=8)
## [1] "Best nuclei c=0.375"
## [1] "95% confidence interval (0.25 0.42)"
#ms <- fitmidparentmodelcs(form="logfitness ~ 1  + distance + (1 | Midparent) + (0+  dummy(Type,\"diploid\")|Acceptor) + (0+  dummy(Type,\"diploid\") |Donor)+ (0+  dummy(Type,\"diploid\")|Acceptor_Donor) + (0+  dummy(Type,\"haploid\") |Donor)", cnucl=seq(0.5, 0.5, by=0.1), type="haploid", relmatrix=list(Acceptor=simmat, Donor=simmat, Midparent=simmat), REML=TRUE, data = data[!(data$Type=="diploid"&data$Acceptor==data$Donor), ])

#ms <- fitmidparentmodelcs(form="logfitness ~ 1  + distance + (1 | Midparent) + (0+  dummy(Type,\"diploid\")|Acceptor) + (0+  dummy(Type,\"diploid\") |Donor)+ (0+  dummy(Type,\"diploid\")|Acceptor_Donor) + (0+  dummy(Type,\"haploid\") |Donor)", cnucl=seq(0.1, 1, by=0.1), type="haploid", relmatrix=list(Acceptor=simmat, Donor=simmat, Midparent=simmat), REML=TRUE, data = data)
summary(ms)
## Linear mixed model fit by REML ['lmerMod']
## 
## REML criterion at convergence: -52212.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.4458 -0.6388 -0.0103  0.6484  3.3599 
## 
## Random effects:
##  Groups         Name                   Variance  Std.Dev. 
##  Acceptor_Donor dummy(Type, "diploid") 2.152e-06 0.0014669
##  Donor          dummy(Type, "haploid") 4.977e-05 0.0070548
##  Donor.1        dummy(Type, "diploid") 3.041e-08 0.0001744
##  Acceptor       dummy(Type, "diploid") 2.741e-08 0.0001656
##  Midparent      (Intercept)            5.729e-03 0.0756919
##  Residual                              1.052e-06 0.0010258
## Number of obs: 4950, groups:  
## Acceptor_Donor, 465; Donor, 30; Acceptor, 30; Midparent, 30
## 
## Fixed effects:
##               Estimate Std. Error t value
## (Intercept) -0.0170554  0.0008735  -19.53
## distance     0.3162111  0.0041707   75.82
## 
## Correlation of Fixed Effects:
##          (Intr)
## distance -0.456

3.3 95%CI

sim <- SimCICompute(simdata=data, nsim=nsimulCI, mbest=ms)
CI <- data.frame(apply(sim, 2, quantile, probs=c(0.025, 0.975)))
write.table(CI, file=paste0('Vars', Sys.Date(),".txt"), sep=" ", row.names = FALSE, quote=NA)

3.4 Check expected and fitted variance among haploid parents

## Check expected and fitted H
Fitted_Var <- data.frame(VarCorr(ms))$vcov[1]

c(Expected_Var=Expected_Var, Fitted_Var=Fitted_Var)
## Expected_Var   Fitted_Var 
## 1.625730e-03 2.151872e-06

3.5 Check expected and fitted H

## Check expected and fitted H
Fitted_Hms=as.numeric(fixef(ms)[2])
c(Expected_H=Expected_Hms, Fitted_H=Fitted_Hms)
## Expected_H   Fitted_H 
##  0.6081143  0.3162111

3.6 Check expected and fitted cnucl

## Check expected and fitted cnucl
Fitted_cnuclms <- attr(ms, "csnucl") 
c(Expected_cnucl=0.5, Fitted_cnucl=Fitted_cnuclms)
## Expected_cnucl   Fitted_cnucl 
##          0.500          0.375
## Compute observed c as ratio of autozygous diploid fitness over haploid fitness 
mean(datadip$meandiplogfitness/datahap$meanhaplogfitness)
## [1] 0.4971465
## Compute observed c as ratio of mean autozygous diploid fitness over mean haploid fitness 
mean(datadip$meandiplogfitness)/mean(datahap$meanhaplogfitness)
## [1] 0.4973536
datadiphapratioms <- datadip$meandiplogfitness/datahap$meanhaplogfitness

## Likelihood profile
#attr(ms, "plotcnucl")

4 Variable h/constant s/constant c: Simulate data with distribution of domiance coefficients across loci (mean h=0.25), but constant s=0.01, and c=0.5

4.1 Prepare data

sim <- simdata(ngeno=ngenosim, nrep=10, Nlocus=Nlocussim, meanshap=0.01, shapesel=NULL, h=0.25, shapeh=2, c=0.5, shapec=NULL, mu=0.0001, sdfit=sdfitsim, seed=seedsim)
## [1] "Simulate data with seed=1"
## [1] "Simulations with distribution of shaploid, h and/or c across selected loci"
## [1] "Mean of dominance coefficients = 0.25"
## [1] "Standard deviation of dominance coefficients = 0.12"
data <- sim$data
distmat <- sim$distmat
Expected_Hmh <- sim$expectedH
datahap <- sim$datahap
datadip <- sim$datadip
Expected_Var <- sim$expectedVar
(observedc <- sim$observedc)
## [1] 0.5
simdat <- 1-distmat
datamh <- data

4.2 Fit model

if(noautozygous){
  data<-data[!(data$Type=="diploid"&data$Acceptor==data$Donor),]
}
#mh <- fitmidparentmodelcs(form="logfitness ~ 1  + distance + (1 | Midparent) + (0+  dummy(Type,\"diploid\")|Acceptor) + (0+  dummy(Type,\"diploid\") |Donor)+ (0+  dummy(Type,\"diploid\")|Acceptor_Donor) + (0+  dummy(Type,\"haploid\") |Donor)", cnucl=seq(0.5, 0.5, by=0.05), type="haploid", relmatrix=NA, REML=TRUE, data = data)

mh <- fitmidparentmodelcs(form="logfitness ~ 1  + distance + (1 | Midparent) + (0+  dummy(Type,\"diploid\")|Acceptor) + (0+  dummy(Type,\"diploid\") |Donor)+ (0+  dummy(Type,\"diploid\")|Acceptor_Donor) + (0+  dummy(Type,\"haploid\") |Donor)", cnucl=seq(start, stop, by=stepc), type="haploid", relmatrix=NA, REML=TRUE, data = data)
## [1] "Fitting model with nuclei cs = 0.1"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26184.03 (df=8)
## [1] "Fitting model with nuclei cs = 0.125"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26186.95 (df=8)
## [1] "Fitting model with nuclei cs = 0.15"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26190.05 (df=8)
## [1] "Fitting model with nuclei cs = 0.175"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26248.98 (df=8)
## [1] "Fitting model with nuclei cs = 0.2"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26173.86 (df=8)
## [1] "Fitting model with nuclei cs = 0.225"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26259.04 (df=8)
## [1] "Fitting model with nuclei cs = 0.25"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26263.38 (df=8)
## [1] "Fitting model with nuclei cs = 0.275"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26267.37 (df=8)
## [1] "Fitting model with nuclei cs = 0.3"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26271.1 (df=8)
## [1] "Fitting model with nuclei cs = 0.325"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26274.68 (df=8)
## [1] "Fitting model with nuclei cs = 0.35"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26278.25 (df=8)
## [1] "Fitting model with nuclei cs = 0.375"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26281.92 (df=8)
## [1] "Fitting model with nuclei cs = 0.4"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26285.84 (df=8)
## [1] "Fitting model with nuclei cs = 0.425"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26290.1 (df=8)
## [1] "Fitting model with nuclei cs = 0.45"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26294.72 (df=8)
## [1] "Fitting model with nuclei cs = 0.475"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26299.47 (df=8)
## [1] "Fitting model with nuclei cs = 0.5"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26303.39 (df=8)
## [1] "Fitting model with nuclei cs = 0.525"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26304.51 (df=8)
## [1] "Fitting model with nuclei cs = 0.55"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26301.46 (df=8)
## [1] "Fitting model with nuclei cs = 0.575"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26202.19 (df=8)
## [1] "Fitting model with nuclei cs = 0.6"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26195.43 (df=8)
## [1] "Fitting model with nuclei cs = 0.625"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26192.27 (df=8)
## [1] "Fitting model with nuclei cs = 0.65"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26187.46 (df=8)
## [1] "Fitting model with nuclei cs = 0.675"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26268.84 (df=8)
## [1] "Fitting model with nuclei cs = 0.7"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26181.56 (df=8)
## [1] "Fitting model with nuclei cs = 0.725"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26179.6 (df=8)
## [1] "Fitting model with nuclei cs = 0.75"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26177.21 (df=8)
## [1] "Fitting model with nuclei cs = 0.775"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26266.26 (df=8)
## [1] "Fitting model with nuclei cs = 0.8"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26266.07 (df=8)
## [1] "Fitting model with nuclei cs = 0.825"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26172.82 (df=8)
## [1] "Fitting model with nuclei cs = 0.85"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26265.74 (df=8)
## [1] "Fitting model with nuclei cs = 0.875"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26265.59 (df=8)
## [1] "Fitting model with nuclei cs = 0.9"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26265.44 (df=8)
## [1] "Fitting model with nuclei cs = 0.925"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26265.31 (df=8)
## [1] "Fitting model with nuclei cs = 0.95"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26265.19 (df=8)
## [1] "Fitting model with nuclei cs = 0.975"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26265.07 (df=8)
## [1] "Fitting model with nuclei cs = 1"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26264.96 (df=8)
## [1] "Best nuclei c=0.525"
## [1] "95% confidence interval (0.44 0.56)"
#mh <- fitmidparentmodelcs(form="logfitness ~ 1  + distance + (1 | Midparent) + (0+  dummy(Type,\"diploid\")|Acceptor) + (0+  dummy(Type,\"diploid\") |Donor)+ (0+  dummy(Type,\"diploid\")|Acceptor_Donor) + (0+  dummy(Type,\"haploid\") |Donor)", cnucl=0.5, type="haploid", relmatrix=list(Acceptor=simmat, Donor=simmat, Midparent=simmat), REML=TRUE, data = data)
summary(mh)
## Linear mixed model fit by REML ['lmerMod']
## 
## REML criterion at convergence: -52609
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3340 -0.6366  0.0009  0.6486  3.3028 
## 
## Random effects:
##  Groups         Name                   Variance  Std.Dev. 
##  Acceptor_Donor dummy(Type, "diploid") 7.811e-07 0.0008838
##  Donor          dummy(Type, "haploid") 2.024e-05 0.0044991
##  Donor.1        dummy(Type, "diploid") 0.000e+00 0.0000000
##  Acceptor       dummy(Type, "diploid") 3.476e-08 0.0001864
##  Midparent      (Intercept)            9.999e-03 0.0999945
##  Residual                              1.052e-06 0.0010258
## Number of obs: 4950, groups:  
## Acceptor_Donor, 465; Donor, 30; Acceptor, 30; Midparent, 30
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept) 0.0053507  0.0009671   5.533
## distance    0.6357843  0.0047321 134.355
## 
## Correlation of Fixed Effects:
##          (Intr)
## distance -0.326

4.3 95%CI

sim <- SimCICompute(simdata=data, nsim=nsimulCI, mbest=mh)
CI <- data.frame(apply(sim, 2, quantile, probs=c(0.025, 0.975)))
write.table(CI, file=paste0('Varh', Sys.Date(),".txt"), sep=" ", row.names = FALSE, quote=NA)

4.4 Check expected and fitted variance among haploid parents

## Check expected and fitted H
Fitted_Var <- data.frame(VarCorr(mh))$vcov[1]
c(Expected_Var=Expected_Var, Fitted_Var=Fitted_Var)
## Expected_Var   Fitted_Var 
## 8.734188e-04 7.811371e-07

4.5 Check expected and fitted H

Fitted_Hmh=as.numeric(fixef(mh)[2])
c(Expected_H=Expected_Hmh, Fitted_H=Fitted_Hmh)
## Expected_H   Fitted_H 
##  0.6178901  0.6357843

4.6 Check expected and fitted cnucl

## Check expected and fitted cnucl
Fitted_cnuclmh <- attr(mh, "csnucl") 
c(Expected_cnucl=0.5, Fitted_cnucl=Fitted_cnuclmh)
## Expected_cnucl   Fitted_cnucl 
##          0.500          0.525
## Compute observed c as ratio of autozygous diploid fitness over haploid fitness 
mean(datadip$meandiplogfitness/datahap$meanhaplogfitness)
## [1] 0.4967349
datadiphapratiomh <- datadip$meandiplogfitness/datahap$meanhaplogfitness

## Likelihood profile
#attr(mh, "plotcnucl")

5 Constant h/constant s/variable c: Simulate data with distribution of ploidy ratios across loci (mean c=0.5), but constant s=0.01, and h=0.25

5.1 Prepare data

sim <- simdata(ngeno=ngenosim, nrep=10, Nlocus=Nlocussim, meanshap=0.01, shapesel=NULL, h=0.25, shapeh=NULL, c=0.5, shapec=2, mu=0.0001, sdfit=sdfitsim, seed=seedsim)
## [1] "Simulate data with seed=1"
## [1] "Simulations with distribution of shaploid, h and/or c across selected loci"
## [1] "Mean of ploidy ratio coefficients = 0.5"
## [1] "Standard deviation of ploidy ratio coefficients = 0.22"
data <- sim$data
distmat <- sim$distmat
Expected_Hmc <- sim$expectedH
datahap <- sim$datahap
datadip <- sim$datadip
Expected_Var <- sim$expectedVar
(observedc <- sim$observedc)
## [1] 0.507092
simdat <- 1-distmat
datamc <- data

5.2 Fit model

if(noautozygous){
  data<-data[!(data$Type=="diploid"&data$Acceptor==data$Donor),]
}
mc <- fitmidparentmodelcs(form="logfitness ~ 1  + distance + (1 | Midparent) + (0+  dummy(Type,\"diploid\")|Acceptor) + (0+  dummy(Type,\"diploid\") |Donor)+ (0+  dummy(Type,\"diploid\")|Acceptor_Donor) + (0+  dummy(Type,\"haploid\") |Donor)", cnucl=seq(start, stop, by=stepc), type="haploid", relmatrix=NA, REML=TRUE, data = data)
## [1] "Fitting model with nuclei cs = 0.1"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26226.26 (df=8)
## [1] "Fitting model with nuclei cs = 0.125"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26229.07 (df=8)
## [1] "Fitting model with nuclei cs = 0.15"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## Warning in optwrap(optimizer, devfun, getStart(start, rho$lower, rho$pp), :
## convergence code 1 from bobyqa: bobyqa -- maximum number of function
## evaluations exceeded
## 'log Lik.' 26232.01 (df=8)
## [1] "Fitting model with nuclei cs = 0.175"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26287.1 (df=8)
## [1] "Fitting model with nuclei cs = 0.2"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26292.21 (df=8)
## [1] "Fitting model with nuclei cs = 0.225"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26296.73 (df=8)
## [1] "Fitting model with nuclei cs = 0.25"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26300.78 (df=8)
## [1] "Fitting model with nuclei cs = 0.275"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26304.44 (df=8)
## [1] "Fitting model with nuclei cs = 0.3"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26307.79 (df=8)
## [1] "Fitting model with nuclei cs = 0.325"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26310.94 (df=8)
## [1] "Fitting model with nuclei cs = 0.35"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26313.99 (df=8)
## [1] "Fitting model with nuclei cs = 0.375"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26316.99 (df=8)
## [1] "Fitting model with nuclei cs = 0.4"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26319.92 (df=8)
## [1] "Fitting model with nuclei cs = 0.425"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26322.68 (df=8)
## [1] "Fitting model with nuclei cs = 0.45"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26324.99 (df=8)
## [1] "Fitting model with nuclei cs = 0.475"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26326.39 (df=8)
## [1] "Fitting model with nuclei cs = 0.5"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26326.37 (df=8)
## [1] "Fitting model with nuclei cs = 0.525"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26324.64 (df=8)
## [1] "Fitting model with nuclei cs = 0.55"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26230.39 (df=8)
## [1] "Fitting model with nuclei cs = 0.575"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26228.41 (df=8)
## [1] "Fitting model with nuclei cs = 0.6"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26223.15 (df=8)
## [1] "Fitting model with nuclei cs = 0.625"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26221.41 (df=8)
## [1] "Fitting model with nuclei cs = 0.65"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26219.65 (df=8)
## [1] "Fitting model with nuclei cs = 0.675"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26217.97 (df=8)
## [1] "Fitting model with nuclei cs = 0.7"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26216.4 (df=8)
## [1] "Fitting model with nuclei cs = 0.725"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26214.95 (df=8)
## [1] "Fitting model with nuclei cs = 0.75"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26214.73 (df=8)
## [1] "Fitting model with nuclei cs = 0.775"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26213.42 (df=8)
## [1] "Fitting model with nuclei cs = 0.8"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26211.35 (df=8)
## [1] "Fitting model with nuclei cs = 0.825"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26210.39 (df=8)
## [1] "Fitting model with nuclei cs = 0.85"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26210.31 (df=8)
## [1] "Fitting model with nuclei cs = 0.875"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26209.51 (df=8)
## [1] "Fitting model with nuclei cs = 0.9"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26300.46 (df=8)
## [1] "Fitting model with nuclei cs = 0.925"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26300.35 (df=8)
## [1] "Fitting model with nuclei cs = 0.95"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26300.25 (df=8)
## [1] "Fitting model with nuclei cs = 0.975"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26300.15 (df=8)
## [1] "Fitting model with nuclei cs = 1"
## [1] "Fitting model with mitochondrial cs = 0"
## [1] "Number of random effects with non-identity variance-covariance matrix:  0"
## 'log Lik.' 26300.06 (df=8)
## [1] "Best nuclei c=0.475"
## [1] "95% confidence interval (0.37 0.5)"
#mc <- fitmidparentmodelcs(form="logfitness ~ 1  + distance + (1 | Midparent) + (0+  dummy(Type,\"diploid\")|Acceptor) + (0+  dummy(Type,\"diploid\") |Donor)+ (0+  dummy(Type,\"diploid\")|Acceptor_Donor) + (0+  dummy(Type,\"haploid\") |Donor)", cnucl=0.5, type="haploid", relmatrix=list(Acceptor=simmat, Donor=simmat, Midparent=simmat), REML=TRUE, data = data)

summary(mc)
## Linear mixed model fit by REML ['lmerMod']
## 
## REML criterion at convergence: -52652.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3517 -0.6357 -0.0003  0.6513  3.2827 
## 
## Random effects:
##  Groups         Name                   Variance  Std.Dev. 
##  Acceptor_Donor dummy(Type, "diploid") 6.750e-07 0.0008216
##  Donor          dummy(Type, "haploid") 5.886e-05 0.0076723
##  Donor.1        dummy(Type, "diploid") 0.000e+00 0.0000000
##  Acceptor       dummy(Type, "diploid") 0.000e+00 0.0000000
##  Midparent      (Intercept)            7.679e-03 0.0876271
##  Residual                              1.052e-06 0.0010258
## Number of obs: 4950, groups:  
## Acceptor_Donor, 465; Donor, 30; Acceptor, 30; Midparent, 30
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept) -0.004623   0.001283  -3.604
## distance     0.625589   0.004590 136.291
## 
## Correlation of Fixed Effects:
##          (Intr)
## distance -0.212

5.3 95%CI

sim <- SimCICompute(simdata=data, nsim=nsimulCI, mbest=mc)
CI <- data.frame(apply(sim, 2, quantile, probs=c(0.025, 0.975)))
write.table(CI, file=paste0('Varc', Sys.Date(),".txt"), sep=" ", row.names = FALSE, quote=NA)

5.4 Check expected and fitted variance among haploid parents

## Check expected and fitted H
Fitted_Var <- data.frame(VarCorr(mc))$vcov[1]
c(Expected_Var=Expected_Var, Fitted_Var=Fitted_Var)
## Expected_Var   Fitted_Var 
## 8.718961e-04 6.749668e-07

5.5 Check expected and fitted H

Fitted_Hmc=as.numeric(fixef(mc)[2])
c(Expected_H=Expected_Hmc, Fitted_H=Fitted_Hmc)
## Expected_H   Fitted_H 
##  0.6338650  0.6255891

5.6 Check expected and fitted cnucl

## Check expected and fitted cnucl
Fitted_cnuclmc <- attr(mc, "csnucl") 
c(Expected_cnuclmc=0.5, Fitted_cnuclmc=Fitted_cnuclmc)
## Expected_cnuclmc   Fitted_cnuclmc 
##            0.500            0.475
## Compute observed c as ratio of autozygous diploid fitness over haploid fitness 
mean(datadip$meandiplogfitness/datahap$meanhaplogfitness)
## [1] 0.4980374
datadiphapratiomc <- datadip$meandiplogfitness/datahap$meanhaplogfitness

## Likelihood profile
#attr(mc, "plotcnucl")

6 Fig. S1: Plot for pH

miny <- min(c(datam00$middipparentheterosis[datam00$Type=="diploid"], datams$middipparentheterosis[datams$Type=="diploid"], datamh$middipparentheterosis[datamh$Type=="diploid"], datamc$middipparentheterosis[datamc$Type=="diploid"]))
maxy <- max(c(datam00$middipparentheterosis[datam00$Type=="diploid"], datams$middipparentheterosis[datams$Type=="diploid"], datamh$middipparentheterosis[datamh$Type=="diploid"], datamc$middipparentheterosis[datamc$Type=="diploid"]))

maxx <- max(c(datam00$distance[datam00$Type=="diploid"], datams$distance[datams$Type=="diploid"], datamh$distance[datamh$Type=="diploid"], datamc$distance[datamc$Type=="diploid"]))
linesp <- 3.5
seglen <- 1

pdf(file="../Figures/FigS1_Fitted_Predicted_H.pdf")

par(mar=c(5.1, 5.1, 4.1, 2.1), mfrow=c(2, 2))
plot(datam00$middipparentheterosis[datam00$Type=="diploid"]~datam00$distance[datam00$Type=="diploid"], xlab="Pairwise genetic distance", ylab=NA, las=1, bty="n", xlim=c(0, maxx), ylim=c(miny, maxy), col="darkgrey")
mtext("Midparent heterosis", side =2, line = linesp)
color <- c("black", "red")
abline(a=0, b=Fitted_Hm00, col=color[1], lwd=2)
#abline(a=0, b=CI$distance[1], col="black", lwd=2, lty=2)
#abline(a=0, b=CI$distance[2], col="black", lwd=2, lty=2)
abline(a=0, b=Expected_Hm00, col=color[2], lwd=2, lty=2)
legend(-0.005, maxy, legend=c("Fitted H", "Predicted H"),lty=c(1, 2), col=color, bty="n", seg.len=seglen, lwd=2)
par(xpd=TRUE)
text(0.01, maxy, "A", cex=2)
par(xpd=FALSE)

plot(datams$middipparentheterosis[datams$Type=="diploid"]~datams$distance[datams$Type=="diploid"], xlab="Pairwise genetic distance", ylab=NA, las=1, bty="n", xlim=c(0, maxx), ylim=c(miny, maxy), col="darkgrey")
mtext("Midparent heterosis", side =2, line = linesp)
abline(a=0, b=Fitted_Hms, col=color[1], lwd=2)
#abline(a=0, b=CI$distance[1], col="black", lwd=2, lty=2)
#abline(a=0, b=CI$distance[2], col="black", lwd=2, lty=2)
abline(a=0, b=Expected_Hms, col=color[2], lwd=2, lty=2)
#legend(0, maxy, legend=c("Predicted H", "Fitted H"),lty=c(1, 2), col=color, bty="n", seg.len=seglen)
par(xpd=TRUE)
text(0.01, maxy, "B", cex=2)
par(xpd=FALSE)


plot(datamh$middipparentheterosis[datamh$Type=="diploid"]~datamh$distance[datamh$Type=="diploid"], xlab="Pairwise genetic distance", ylab=NA, las=1, bty="n", xlim=c(0, maxx), ylim=c(miny, maxy), col="darkgrey")
mtext("Midparent heterosis", side =2, line = linesp)
abline(a=0, b=Fitted_Hmh, col=color[1], lwd=2)
#abline(a=0, b=CI$distance[1], col="black", lwd=2, lty=2)
#abline(a=0, b=CI$distance[2], col="black", lwd=2, lty=2)
abline(a=0, b=Expected_Hmh, col=color[2], lwd=2, lty=2)
#legend(0, maxy, legend=c("Predicted H", "Fitted H"),lty=c(1, 2), col=color, bty="n", seg.len=seglen)
par(xpd=TRUE)
text(0.01, maxy, "C", cex=2)
par(xpd=FALSE)


plot(datamc$middipparentheterosis[datamc$Type=="diploid"]~datamc$distance[datamc$Type=="diploid"], xlab="Pairwise genetic distance", ylab=NA, las=1, bty="n", xlim=c(0, maxx), ylim=c(miny, maxy), col="darkgrey")
mtext("Midparent heterosis", side =2, line = linesp)
abline(a=0, b=Fitted_Hmc, col=color[1], lwd=2)
#abline(a=0, b=CI$distance[1], col="black", lwd=2, lty=2)
#abline(a=0, b=CI$distance[2], col="black", lwd=2, lty=2)
abline(a=0, b=Expected_Hmc, col=color[2], lwd=2, lty=2)
#legend(0, maxy, legend=c("Predicted H", "Fitted H"), lty=c(1, 2), col=color, bty="n", seg.len=seglen)
par(xpd=TRUE)
text(0.01, maxy, "D", cex=2)


dev.off()
## quartz_off_screen 
##                 2

7 Fig.S2: Histogram for cnucl

maxy <- 0.8
lety <- maxy 
#brks <- seq(round(min(range(c(datam00, datams, datamh, datamc))), 0)-1, round(max(range(c(datam00, datams, datamh, datamc))), 0)+1, by = 1)

minx <- -0.5
maxx <- 1.5

letx <- -0.4
seglen <- 1
datadiphapratiom00 <- datadiphapratiom00[datadiphapratiom00>minx&datadiphapratiom00<maxx]
datadiphapratioms <- datadiphapratioms[datadiphapratioms>minx&datadiphapratioms<maxx]
datadiphapratiomh <- datadiphapratiomh[datadiphapratiomh>minx&datadiphapratiomh<maxx]
datadiphapratiomc <- datadiphapratiomc[datadiphapratiomc>minx&datadiphapratiomc<maxx]

brks <- seq(minx, maxx, by=0.05)
pdf(file="../Figures/FigS2_Expected_Fitted_cnucl.pdf")

par(mfrow=c(2, 2))
#par(mar=c(5, 5, 2, 1), mfrow=c(2, 2), xpd=TRUE)

histFull <- hist(datadiphapratiom00, breaks = brks, plot = FALSE)
histFull$counts <- histFull$counts/sum(histFull$counts)
plot(histFull, col = "darkgrey", xlab = "Fitness ratio of autozygous diploids over haploids", ylab = "Frequency distribution", axes = TRUE, main = "", freq=TRUE, xlim=range(brks), ylim=c(0, maxy), las=1, border = "white")
abline(v=0.5, lty=2, col="red")
abline(v=Fitted_cnuclm00, lty=1, col="black")
color <- c("black", "red")

legend(-0.6, maxy, legend=c("Fitted H", "Predicted H"),lty=c(1, 2), col=color, bty="n", seg.len=seglen, lwd=2)

par(xpd=TRUE)
text(letx, lety, "A", cex=2)
par(xpd=FALSE)


histFull <- hist(datadiphapratioms, breaks = brks, plot = FALSE)
histFull$counts <- histFull$counts/sum(histFull$counts)
plot(histFull, col = "darkgrey", xlab = "Fitness ratio of autozygous diploids over haploids", ylab = "Frequency distribution", axes = TRUE, main = "", freq=TRUE, xlim=range(brks), ylim=c(0, maxy), las=1, border = "white")
abline(v=0.5, lty=2, col="red")
abline(v=Fitted_cnuclms, lty=1, col="black")

par(xpd=TRUE)
text(letx, lety, "B", cex=2)
par(xpd=FALSE)

histFull <- hist(datadiphapratiomh, breaks = brks, plot = FALSE)
histFull$counts <- histFull$counts/sum(histFull$counts)
plot(histFull, col = "darkgrey", xlab = "Fitness ratio of autozygous diploids over haploids", ylab = "Frequency distribution", axes = TRUE, main = "", freq=TRUE, xlim=range(brks), ylim=c(0, maxy), las=1, border = "white")
abline(v=0.5, lty=2, col="red")
abline(v=Fitted_cnuclmh, lty=1, col="black")

par(xpd=TRUE)
text(letx, lety, "C", cex=2)
par(xpd=FALSE)

histFull <- hist(datadiphapratiomc, breaks = brks, plot = FALSE)
histFull$counts <- histFull$counts/sum(histFull$counts)
plot(histFull, col = "darkgrey", xlab = "Fitness ratio of autozygous diploids over haploids", ylab = "Frequency distribution", axes = TRUE, main = "", freq=TRUE, xlim=range(brks), ylim=c(0, maxy), las=1, border = "white")
abline(v=0.5, lty=2, col="red")
abline(v=Fitted_cnuclmc, lty=1, col="black")

par(xpd=TRUE)
text(letx, lety, "D", cex=2)
par(xpd=FALSE)

dev.off()
## quartz_off_screen 
##                 2