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
Constant h/s/c: Simulate data with h=0.25, s=0.01 and c=0.5
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
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
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)
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
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
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")
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
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
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
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)
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
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
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")
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
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
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
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)
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
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
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")
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
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
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
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)
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
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
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")
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
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