version 0.1.1

cran · Sep 6, 2016 · d627774 · d627774
commit d627774
Show file tree

Hide file tree

Showing 96 changed files with 10,910 additions and 0 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -0,0 +1,16 @@
+Package: SpadeR
+Type: Package
+Title: Species-Richness Prediction and Diversity Estimation with R
+Version: 0.1.1
+Date: 2016-09-06
+Author: Anne Chao, K. H. Ma, T. C. Hsieh and Chun-Huo Chiu
+Maintainer: Anne Chao <chao@stat.nthu.edu.tw>
+Description: Estimation of various biodiversity indices and related (dis)similarity measures based on individual-based (abundance) data or sampling-unit-based (incidence) data taken from one or multiple communities/assemblages.
+License: GPL (>= 3)
+Imports: stats
+Depends: R (>= 2.14)
+RoxygenNote: 5.0.1
+NeedsCompilation: no
+Packaged: 2016-09-06 05:40:30 UTC; USER
+Repository: CRAN
+Date/Publication: 2016-09-06 14:53:57
diff --git a/MD5 b/MD5
@@ -0,0 +1,95 @@
+ef63b6781a5d00b24a8ceab08ec09fe7 *DESCRIPTION
+cc5423c19ec4d8c4dbe8294ed31230ef *NAMESPACE
+96457fe9e2885e62a1730637e7585d58 *R/BasicFun.R
+5d7753da40d64deb34cafe045c765a23 *R/BootstrapFunMa.R
+c334bc7b849cf6539ecb451226c9df20 *R/Cf0Fun.R
+cb0265a58940b0cce40a390f757bcbd0 *R/Chao1_bcEstFun.R
+b3fb7c9989c503646096f1abd50a6547 *R/Chao1_bcFun.R
+039be360ad01867f44254330949c6e1c *R/Chao1_sharedEstFun.R
+7c4a548deef925116e8a48984a1addc1 *R/Chao1_sharedFun.R
+3f77a1cdcf2f1b2fb5263817f44a8f88 *R/Chao2_bcEstFun.R
+c153837bbc449206f7ea681e6d94fd28 *R/Chao2_bcFun.R
+c89b36add1fbc292d20f6eb55c2ce85f *R/Chao2_sharedEstFun.R
+cdd0cbf8f385812ff5f55794227e5a5d *R/Chao2_sharedFun.R
+d7469eabd5f4cf5b483c553ad6a8fe4e *R/ChaoShared.Ind.R
+ca71c23b7a874f01ffdf1a6dd33c3696 *R/ChaoShared.Sam.R
+cbee1eb899e2ddc829cc2ef62fe93748 *R/DataTransform.R
+fe89e0e7eae2864010f4f9a7f6f65909 *R/Diversity_subroutine.R
+6e5bd623f03da44acdf0d4540fa18cd0 *R/ExtenProbFunMa.R
+d980172e8d892a76caf5a2cf79f19335 *R/Genetic_subroutine.R
+c00bbe6954c76b4b27602fe3bc288439 *R/HeteroEstFun.R
+f1f0e4cf7ecb19fd7412805156ed8446 *R/HeteroFun.R
+83517ffa83e8b8e087904f10c9168b78 *R/HomoEstFun.R
+3d9a45bb13ce034198be5d1cf71c86fb *R/HomoFun.R
+d8dfae8845816236c93c32cdcb813438 *R/InfreqSpeciesGroupInprove.R
+987f73637aa92efcfaf395505d6ccf82 *R/Multiple_Community_Measure_subroutine.R
+39dc8325397490ed3140fe8a0eb115e7 *R/PanEstFun.R
+89a11bdd038507caf5f60c5485125f9b *R/PanEstFun.Sam.R
+4769699365e75ef21c150f2240d06ca6 *R/PanFun.R
+50c6687b94e782758c0fa10dc05b165a *R/PanFun.Sam.R
+d18b8fe73323d4a9c371a5bf8141c777 *R/PanbcEstFun.R
+a4a5d47abe22d13a1e5c38a0864a5949 *R/PanbcEstFun.Sam.R
+833b94bcf514a054b073792ad26acd28 *R/PanbcFun.R
+d90ec73de3b2b91288b73276754f85dd *R/PanbcFun.Sam.R
+1b0cf8a5e1f11a53d261294b3dc577ff *R/Q.R
+cd3ddba9f74ce2de02a929414d52e85b *R/RareSpeciesGroup.R
+337792c0859f55b0920aa893ae7e785b *R/RareSpeciesGroupInprove.R
+a5682eb9cad07cef9939ea16619c6f9c *R/SortDataFun.R
+5bd3709d9e635e97d18bb61cdb116680 *R/SpecAbunAce.R
+79e74c922e2c0598fc10aefba21974fc *R/SpecAbunAce1.R
+7e41eb21e99a160367481f9744c8cd60 *R/SpecAbunChao1.R
+ab21bf01562f267d0f7d11c6f02719fd *R/SpecAbunChao1bc.R
+f62cbc6946fe98b5e36830a0844b3078 *R/SpecAbunHomo.R
+b21e53f1c4ca1840c5103d0f18dfa2ac *R/SpecAbunHomoMle.R
+42390b87ed8cae0a0f8db06c69b1224d *R/SpecAbunJack1.R
+681da6a5f421d9013dd1e75b958bb871 *R/SpecAbunJack2.R
+30ccbe67d4577795408bbd8bac8a84a2 *R/SpecAbunOut.R
+a51e167166a6d1d17b8bb3d0542942c5 *R/SpecAbuniChao1.R
+c7227301d3b6fa7a5b115354f380b142 *R/SpecInciChao2.R
+803be732f3dcf3100349318b899a2c33 *R/SpecInciChao2bc.R
+b669edf7d22082ec617e7271384dc86f *R/SpecInciHomo.R
+45b25c518af62ce91a89f2762a725ea0 *R/SpecInciJack1.R
+87f660e1cbee3bfbd7d94ab90de835cb *R/SpecInciJack2.R
+867c837cf023c8003bc314b3227ab7da *R/SpecInciModelh.R
+593fcd8e64681d9737d4d9106a13ae1a *R/SpecInciModelh1.R
+5e43e5933dea74aacad4469b1d0f7649 *R/SpecInciModelth.R
+aa944c0a60bf41c64aaf797f72ca579b *R/SpecInciModelth1.R
+482486ae9d1f4e36466e4f005f78f0e5 *R/SpecInciOut.R
+bd8f9a1282d1e3df2048307f13057da6 *R/SpecInciiChao2.R
+4ca5e2eeabaac55241e63b572870c3a7 *R/SpeciesAbundance.R
+a320dc46f58c26278c4ef7cf6d6f6fb2 *R/SpeciesIncidence.R
+cf42cbba7113351d64dd37f8dbbdbf52 *R/Two_Community_similarity_subroutine.R
+c0d0e1cc1cb39f11182ad8520967f128 *R/VarEstFun.R
+e649a92435617a44bceabb8516a522ed *R/VarEstFun.Sam.R
+2a0b9fa3376a386b4c69da5668aacb44 *R/basicAbun.R
+f593c24ad2d29cc948107730f5c3770d *R/basicInci.R
+1c1ff82971d608d0fe05d66563a56483 *R/diff_Chao1.R
+8b08b749f3ee692d1f260f9a736f2d35 *R/diff_Chao1bc.R
+8b324719c0b5d1e10ae7ca78224224fd *R/diff_Chao2.R
+68e1fd85b96535239af6c0516cd98357 *R/diff_Chao2bc.R
+908c46104df13696e05df6b165e4fdf5 *R/diff_Pan.R
+7452f504e68c93f2bfd9eae8f579bdeb *R/diff_Panbc.R
+218fd6e8ce56b30d220b8bdc3bf6d315 *R/f.R
+a5744186fadffea992bfc37f66b43974 *R/logCI.R
+f13660ab18e0d0274707d659badfbd71 *R/print.ChaoShared.R
+be07b5b687c48cc73df8c5f3ccd3debd *R/print.ChaoSpecies.R
+d77cef44b1a076fa6719fcda59ad422a *R/spader.R
+0ac3c5357791fc3998def99e6c0f7c82 *data/ChaoSharedData.rda
+b78e7681dd20cd52f35879ba61a6cb4f *data/ChaoSpeciesData.rda
+52ca0aaa3349ee6db5ec9e241567a5fd *data/DiversityData.rda
+0352f4e62caada865d1cbcf9f5ed27e4 *data/GeneticsDataAbu.rda
+5a832c3b55c6048eecf020e50a0a591a *data/SimilarityMultData.rda
+b938ca1206c477d439364ea5ac20a227 *data/SimilarityPairData.rda
+b911af0e3d807849a3d4ac2143bb216d *man/ChaoShared.Rd
+842d4a37df37aa2ebff1507f2a811f46 *man/ChaoSpecies.Rd
+2d08ebb6103373a38ed376afa3f4ecba *man/DivDemoData.Rd
+8c3205b588f99e61d2b8e15fb5af1988 *man/Diversity.Rd
+960f74e6b04717faa5f1f2ac331fcde6 *man/Genetics.Rd
+ef85ac24c7e724ab0c7caea534c9c486 *man/GeneticsDataAbu.Rd
+b0fa6c87e38b26e832ccb1b8c7d12ab3 *man/MulComSimDemoData.Rd
+c5af1479718c074e02cd3b98d756956d *man/SharedSpecDemoData.Rd
+5041135dca36d2b5cb4f413c60694ad3 *man/SimilarityMult.Rd
+b685a0965a551f5775f8360e92e258a6 *man/SimilarityPair.Rd
+921b76b91b118770fc93c3cf1641f408 *man/SpecDemoData.Rd
+a1584e2d81563c4f968d707df9eced61 *man/TwoComSimDemoData.Rd
+96c41d4f4d6b4fe6866ca6893674cf7a *man/spader-package.Rd
diff --git a/NAMESPACE b/NAMESPACE
@@ -0,0 +1,15 @@
+import(stats)
+importFrom("grDevices", "adjustcolor")
+importFrom("graphics", "legend", "lines", "plot", "polygon")
+S3method(print,ChaoShared)
+S3method(print,ChaoSpecies)
+S3method(print,spadeDiv)
+S3method(print,spadeTwo)
+S3method(print,spadeMult)
+S3method(print,spadeGenetic)
+export(ChaoShared)
+export(ChaoSpecies)
+export(Diversity)
+export(SimilarityMult)
+export(SimilarityPair)
+export(Genetics)
diff --git a/R/BasicFun.R b/R/BasicFun.R
@@ -0,0 +1,81 @@
+BasicFun <-  function(x1, x2, B, datatype) {
+  if(datatype=="abundance"){
+    n1 <- sum(x1)
+    n2 <- sum(x2)
+    D1 <- sum(x1 > 0)
+    D2 <- sum(x2 > 0)
+    D12 <- sum(x1 > 0 & x2 > 0)
+    x1_share <- x1[which(x1 > 0 & x2 > 0)]
+    x2_share <- x2[which(x1 > 0 & x2 > 0)]
+    f11 <- sum(x1_share == 1 & x2_share == 1)
+    f1.plus.rare <- sum(x1_share == 1 & x2_share <= 10)
+    fplus.1.rare <- sum(x2_share == 1 & x1_share <= 10)
+    f2.plus.rare <- sum(x1_share == 2 & x2_share <= 10)
+    fplus.2.rare <- sum(x2_share == 2 & x1_share <= 10)
+    f1.plus <- sum(x1_share == 1 & x2_share >0)
+    fplus.1 <- sum(x2_share == 1 & x1_share >0)
+    f2.plus <- sum(x1_share == 2 & x2_share >0)
+    fplus.2 <- sum(x2_share == 2 & x1_share >0)
+    f22 <- sum(x1_share == 2 & x2_share == 2)
+    D12_rare <- sum(x1_share <= 10 & x2_share <= 10)
+
+    pos <- (x1 > 0 & x2 > 0) & (x1 > 10 | x2 > 10)
+    n1_rare <- n1 - sum(x1[pos])
+    n2_rare <- n2 - sum(x2[pos])
+
+    pos_r <- (x1_share <= 10 & x2_share <= 10)
+    pos1_r <- (x1_share == 1 & x2_share <= 10)
+    pos2_r <- (x2_share == 1 & x1_share <= 10)
+
+    tmp <- sum(x1_share[pos_r] * x2_share[pos_r])
+    C12_rare <- 1 - (sum(x2_share[pos1_r]) + sum(x1_share[pos2_r]) - f11) / tmp
+    #   C12_rare <- round(C12_rare, 4)
+
+    T10 <- sum(x1_share[x1_share <= 10 & x2_share <= 10])
+    T01 <- sum(x2_share[x1_share <= 10 & x2_share <= 10])
+    T11 <- tmp
+    T21 <- sum(x1_share[pos_r] * (x1_share - 1)[pos_r] * x2_share[pos_r])
+    T12 <- sum(x1_share[pos_r] * (x2_share - 1)[pos_r] * x2_share[pos_r])
+
+    T22 <- sum(x1_share[pos_r] * x2_share[pos_r] * 
+                 (x1_share - 1)[pos_r] * (x2_share - 1)[pos_r])
+
+    S12_0 <- D12_rare / C12_rare
+    CCV_1 <- S12_0 * n1_rare * T21 / (n1_rare - 1) / T10 / T11 - 1
+    CCV_2 <- S12_0 * n2_rare * T12 / (n2_rare - 1) / T01 / T11 - 1
+    CCV_12 <- n1_rare * n2_rare * S12_0^2 * T22 / 
+      ((n1_rare - 1) * (n2_rare - 1) * T10 * T01 * T11) - 
+      S12_0 * T11 / T10 / T01 - CCV_1 - CCV_2
+    out <- list(n1=n1, n2=n2, D1=D1, D2=D2, D12=D12, B=B, f11=f11, f1.plus=f1.plus,
+                 fplus.1=fplus.1, f2.plus=f2.plus, fplus.2=fplus.2, f22=f22,
+                f1.plus.rare=f1.plus.rare,fplus.1.rare=fplus.1.rare, f2.plus.rare=f2.plus.rare, fplus.2.rare=fplus.2.rare,
+                n1_rare=n1_rare, n2_rare=n2_rare, D12_rare=D12_rare, 
+                C12_rare=C12_rare, CCV_1=CCV_1, CCV_2=CCV_2, CCV_12=CCV_12,
+                datatype="abundance")
+
+    }else if(datatype=="incidence_freq"){
+    y1 <- x1
+    y2 <- x2
+    t1 <- y1[1]
+    t2 <- y2[1]
+    x1 <- y1[-1]
+    x2 <- y2[-1]
+    u1 <- sum(x1)
+    u2 <- sum(x2)
+    D1 <- sum(x1 > 0)
+    D2 <- sum(x2 > 0)
+    D12 <- sum(x1 > 0 & x2 > 0)  
+
+    Q11 <- sum(x1 == 1 & x2 == 1)
+    Q1.plus <- sum(x1 == 1 & x2 >= 1)
+    Qplus.1 <- sum(x2 == 1 & x1 >= 1)
+    Q2.plus <- sum(x1 == 2 & x2 >= 1)
+    Qplus.2 <- sum(x2 == 2 & x1 >= 1)
+    Q22 <- sum(x1 == 2 & x2 == 2) 
+    out <- list(T1=t1, T2=t2, U1=u1, U2=u2, D1=D1, D2=D2, D12=D12, B=B, Q11=Q11, Q1.plus=Q1.plus,
+                 Qplus.1=Qplus.1, Q2.plus=Q2.plus, Qplus.2=Qplus.2, Q22=Q22,
+                datatype="incidence")  
+
+    }else warning("invaild datatype, please set datatype as abundance or incidence.")
+    return(out)
+  }
diff --git a/R/BootstrapFunMa.R b/R/BootstrapFunMa.R
@@ -0,0 +1,26 @@
+BootstrapFunMa <-
+function(x1, x2, B, FunName) {
+  n1 <- sum(x1); n2 <- sum(x2)
+  z <- SortDataFun(x1, x2)
+  z1 <- z[, 1]
+  z2 <- z[, 2]
+  newprob <- ExtenProbFunMa(z1, z2)
+  p1 <- newprob$prob1
+  p2 <- newprob$prob2
+  set.seed(123)
+  X1 <- rmultinom(B, n1, p1)
+  set.seed(123)
+  X2 <- rmultinom(B, n2, p2)
+  X <- rbind(X1, X2)
+
+  se <- sd(apply(X, 2, function(x) {
+    y <- matrix(x, ncol=2)
+    y1 <- y[, 1]
+    y2 <- y[, 2]
+#     y1 <- x[1 : length(p1)]
+#     y2 <- x[(length(p1) + 1) : (2 * length(p1))]
+    FunName(y1, y2)
+  }), na.rm=T)
+
+  return(se)
+}
diff --git a/R/Cf0Fun.R b/R/Cf0Fun.R
@@ -0,0 +1,17 @@
+Cf0Fun <-
+function(x) {
+  n <- sum(x)
+  f1 <- sum(x == 1); f2 <- sum(x == 2)
+  if (f2 > 0) {
+    C <- 1 - f1 / n * ((n - 1) * f1 / ((n - 1) * f1 + 2 * f2))
+    f0 <- (n - 1) / n * f1^2 / (2 * f2)
+  } else if (f2 == 0 & f1 != 0) {
+    C <- 1 - f1 / n * ((n - 1) * (f1 - 1) / ((n - 1) * (f1 - 1) + 2))
+    f0 <- (n - 1) / n * f1 * (f1 - 1) / 2
+  } else {
+    C <- 1
+    f0 <- (n - 1) / n * f1 * (f1 - 1) / 2
+  }
+  f0 <- ceiling(f0)
+  return(c(C, f0))
+}
diff --git a/R/Chao1_bcEstFun.R b/R/Chao1_bcEstFun.R
@@ -0,0 +1,15 @@
+Chao1_bcEstFun <-
+function(x1, x2) {
+  D12 <- sum(x1 > 0 & x2 > 0)
+  x1_share <- x1[which(x1 > 0 & x2 > 0)]
+  x2_share <- x2[which(x1 > 0 & x2 > 0)]
+  f11 <- sum(x1_share == 1 & x2_share == 1)
+  f1.plus <- sum(x1_share == 1 & x2_share >= 1)
+  fplus.1 <- sum(x2_share == 1 & x1_share >= 1)
+  f2.plus <- sum(x1_share == 2 & x2_share >= 1)
+  fplus.2 <- sum(x2_share == 2 & x1_share >= 1)
+  est <- D12 + f11 * f1.plus * fplus.1 / (4 * (f2.plus + 1) * (fplus.2 + 1)) + 
+    f1.plus * (f1.plus - 1) / (2 * (f2.plus + 1)) + 
+    fplus.1 * (fplus.1 - 1) / (2 * (fplus.2 + 1))
+  return(est)
+}
diff --git a/R/Chao1_bcFun.R b/R/Chao1_bcFun.R
@@ -0,0 +1,11 @@
+Chao1_bcFun <-
+function(x1, x2, conf=0.95) {
+  est <- Chao1_bcEstFun(x1, x2)
+  se <- VarEstFun(x1, x2, diffFun=diff_Chao1bc, FunName=Chao1_bcEstFun)
+  CI <- logCI(x1, x2, est, se, conf)
+  out <- matrix(c(est, se, CI), nrow = 1)
+  rownames(out) <- c("Chao1-shared-bc")
+  colnames(out) <- c("Estimator", "Est_s.e.",
+                     paste(conf*100, "% Lower"), paste(conf*100, "% Upper"))
+  return(out)
+}
diff --git a/R/Chao1_sharedEstFun.R b/R/Chao1_sharedEstFun.R
@@ -0,0 +1,20 @@
+Chao1_sharedEstFun <-
+function(x1, x2) {
+  D12 <- sum(x1 > 0 & x2 > 0)
+  x1_share <- x1[which(x1 > 0 & x2 > 0)]
+  x2_share <- x2[which(x1 > 0 & x2 > 0)]
+  f11 <- sum(x1_share == 1 & x2_share == 1)
+  f1.plus <- sum(x1_share == 1 & x2_share >= 1)
+  fplus.1 <- sum(x2_share == 1 & x1_share >= 1)
+  f2.plus <- sum(x1_share == 2 & x2_share >= 1)
+  fplus.2 <- sum(x2_share == 2 & x1_share >= 1)
+  if (f2.plus == 0 || fplus.2 == 0) {
+    est <- D12 + f11 * f1.plus * fplus.1 / (4 * (f2.plus + 1) * (fplus.2 + 1)) + 
+      f1.plus * (f1.plus - 1) / (2 * (f2.plus + 1)) + 
+      fplus.1 * (fplus.1 - 1) / (2 * (fplus.2 + 1))
+  } else {
+    est <- D12 + f11 * f1.plus * fplus.1 / (4 * f2.plus * fplus.2) + 
+      f1.plus^2 / (2 * f2.plus) + fplus.1^2 / (2 * fplus.2)
+  }
+  return(est)
+}
diff --git a/R/Chao1_sharedFun.R b/R/Chao1_sharedFun.R
@@ -0,0 +1,19 @@
+Chao1_sharedFun <-
+function(x1, x2, conf=0.95) {
+  f2p <- sum(x1 == 2 & x2 >= 1)
+  fp2 <- sum(x1 >= 1 & x2 == 2)
+
+  est <- Chao1_sharedEstFun(x1, x2)
+  if (f2p == 0 || fp2 == 0) { 
+    se <- VarEstFun(x1, x2, diffFun=diff_Chao1bc, FunName=Chao1_bcEstFun)
+  } else {
+    se <- VarEstFun(x1, x2, diff_Chao1, FunName=Chao1_sharedEstFun)
+  }
+
+  CI <- logCI(x1, x2, est, se, conf)
+  out <- matrix(c(est, se, CI), nrow = 1)
+  rownames(out) <- c("Chao1(shared)")
+  colnames(out) <- c("Estimator", "Est_s.e.",
+                     paste(conf*100, "% Lower"), paste(conf*100, "% Upper"))
+  return(out)
+}
diff --git a/R/Chao2_bcEstFun.R b/R/Chao2_bcEstFun.R
@@ -0,0 +1,19 @@
+Chao2_bcEstFun <-
+function(y1, y2) {
+  t1 <- y1[1]
+  t2 <- y2[1]
+  c1 <- (t1 - 1) / t1
+  c2 <- (t2 - 1) / t2
+  x1 <- y1[-1]
+  x2 <- y2[-1]
+  D12 <- sum(x1 > 0 & x2 > 0)
+  Q11 <- sum(x1 == 1 & x2 == 1)
+  Q1.plus <- sum(x1 == 1 & x2 >= 1)
+  Qplus.1 <- sum(x2 == 1 & x1 >= 1)
+  Q2.plus <- sum(x1 == 2 & x2 >= 1)
+  Qplus.2 <- sum(x2 == 2 & x1 >= 1)
+  est <- D12 + Q11 * c1 * c2 * Q1.plus * Qplus.1 / (4 * (Q2.plus + 1) * (Qplus.2 + 1)) + 
+    c1 * Q1.plus * (Q1.plus - 1) / (2 * (Q2.plus + 1)) + 
+    c2 * Qplus.1 * (Qplus.1 - 1) / (2 * (Qplus.2 + 1))
+  return(est)
+}
diff --git a/R/Chao2_bcFun.R b/R/Chao2_bcFun.R
@@ -0,0 +1,11 @@
+Chao2_bcFun <-
+function(y1, y2, conf=0.95) {
+  est <- Chao2_bcEstFun(y1, y2)
+  se <- VarEstFun.Sam(y1, y2, diffFun=diff_Chao2bc, FunName=Chao2_bcEstFun)
+  CI <- logCI(y1, y2, est, se, conf)
+  out <- matrix(c(est, se, CI), nrow = 1)
+  rownames(out) <- c("Chao2-shared-bc")
+  colnames(out) <- c("Estimator", "Est_s.e.",
+                     paste(conf*100, "% Lower"), paste(conf*100, "% Upper"))
+  return(out)
+}
diff --git a/R/Chao2_sharedEstFun.R b/R/Chao2_sharedEstFun.R
@@ -0,0 +1,24 @@
+Chao2_sharedEstFun <-
+function(y1, y2) {
+  t1 <- y1[1]
+  t2 <- y2[1]
+  c1 <- (t1 - 1) / t1
+  c2 <- (t2 - 1) / t2
+  x1 <- y1[-1]
+  x2 <- y2[-1]
+  D12 <- sum(x1 > 0 & x2 > 0)  
+  Q11 <- sum(x1 == 1 & x2 == 1)
+  Q1.plus <- sum(x1 == 1 & x2 >= 1)
+  Qplus.1 <- sum(x2 == 1 & x1 >= 1)
+  Q2.plus <- sum(x1 == 2 & x2 >= 1)
+  Qplus.2 <- sum(x2 == 2 & x1 >= 1)  
+  if (Q2.plus == 0 || Qplus.2 == 0) {
+    est <- D12 + Q11 * c1 * c2 * Q1.plus * Qplus.1 / (4 * (Q2.plus + 1) * (Qplus.2 + 1)) + 
+      c1 * Q1.plus * (Q1.plus - 1) / (2 * (Q2.plus + 1)) + 
+      c2 * Qplus.1 * (Qplus.1 - 1) / (2 * (Qplus.2 + 1))
+  } else {
+    est <- D12 + Q11 * c1 * c2 * Q1.plus * Qplus.1 / (4 * Q2.plus * Qplus.2) + 
+      c1 * Q1.plus^2 / (2 * Q2.plus) + c2 * Qplus.1^2 / (2 * Qplus.2)
+  }
+  return(est)
+}