version 1.0.4

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: likelihoodR
 Title: Likelihood Analyses for Common Statistical Tests
-Version: 1.0.3
+Version: 1.0.4
 Authors@R: 
     person(given = "Peter",
            family = "Cahusac",
@@ -20,11 +20,11 @@ Description: A collection of functions that calculate the log likelihood
     P. Cahusac (2020, ISBN-13 : 978-1119549802).
 License: GPL-2
 Encoding: UTF-8
-LazyData: true
+LazyData: false
 RoxygenNote: 7.1.1
 NeedsCompilation: no
-Packaged: 2021-02-07 11:17:49 UTC; pcahusac
+Packaged: 2021-05-03 12:16:01 UTC; pcahusac
 Author: Peter Cahusac [aut, cre] (<https://orcid.org/0000-0003-4976-2834>)
 Maintainer: Peter Cahusac <peteqsac@gmail.com>
 Repository: CRAN
-Date/Publication: 2021-02-11 09:50:02 UTC
+Date/Publication: 2021-05-03 12:30:02 UTC

MD5

@@ -1,26 +1,26 @@
-65532fe19de0da1b5945cd3a8e5ac728 *DESCRIPTION
-23f6d5270a3fd68dd48f41563282f3c4 *NAMESPACE
-7b3eb24ddca0b0eeeb3324a4283f8924 *R/L_1way_ANOVA.R
-0d0eac744ada49d0b3ea8805f0adaa5e *R/L_1way_RM_ANOVA.R
-84d0534bd58fd3e93f5ba996a325f7f4 *R/L_1way_cat.R
+ce114942a9f9a3ae12f8ba70f2c1195c *DESCRIPTION
+a17a73059d9ce8de56237528a0b3c965 *NAMESPACE
+5c55bd2e51514e28c91a13acc2faf1fb *R/L_1way_ANOVA.R
+5c03a9613ee4fb2e7ef1e321b8c2a599 *R/L_1way_RM_ANOVA.R
+0262779a5a0822be8d772eed6251e38c *R/L_1way_cat.R
 22aaf99cf48f11110412a72c88e7f456 *R/L_2S_ttest.R
-f1d2170b479fd4cdc32d78182a89f8e1 *R/L_2way_Factorial_ANOVA.R
-1e7f55bf50b033f076041c99e6b4efa5 *R/L_2way_cat.R
-1cce946a0a70230980958abfc2563317 *R/L_OR.R
-138b9eb1ae6abdbe03acfd65e1e8f6a3 *R/L_RR.R
+5bb61c75c922e04646309a490ae43aa5 *R/L_2way_Factorial_ANOVA.R
+161048edea31e3721d5f8acd5b9074ef *R/L_2way_cat.R
+3e5aa97fe80c8055a21c73b0e11c8243 *R/L_OR.R
+ba3231359ae369c75f2cb6fe61d0e9dd *R/L_RR.R
 9dfa30e0b7c8cc7ad2f4945b751b0c94 *R/L_corr.R
-d6430f86082570e6329b8cc5347366ea *R/L_efficacy.R
+cf8fbfa2301deb25dcc7bf30bbeef942 *R/L_efficacy.R
 b5bd708114dabb8138c7273a563c83d3 *R/L_logistic_regress.R
-7aace07b72e173e2f17f1c23abe71cae *R/L_regress.R
+89199e2bb5cb0982b16d83c574da5c8e *R/L_regress.R
 dd4f02f6e6a9f896a06cfa9bf604f8b3 *R/L_t_test_sample_size.R
 65c1853af5f1489bd3ee292fa42c8d02 *R/L_ttest.R
 20350b22f3ccda17c9795d70d93e71aa *man/L_1way_ANOVA.Rd
 145baf4e6ffe0a65c3fddb5706244f13 *man/L_1way_RM_ANOVA.Rd
-4575a09aa8fc4a7bfd7753c5fea50835 *man/L_1way_cat.Rd
+d61c424bef4624993ee48d1a1bdbf69e *man/L_1way_cat.Rd
 2b483ca9b7085c8563e820fb90890535 *man/L_2S_ttest.Rd
-f744d69e0ce7ef58dd23682196362268 *man/L_2way_Factorial_ANOVA.Rd
-9217e4f57e3add4ac9dbe3b230959911 *man/L_2way_cat.Rd
-be180e5d0ba382039ba0c66511bd6b82 *man/L_OR.Rd
+d15cb8dec9a4916daa594c6432538fd5 *man/L_2way_Factorial_ANOVA.Rd
+a252bd17c6d53fe14cf9796e88f0d3a1 *man/L_2way_cat.Rd
+9d56afc1bd44707c4c3c73686fe7642c *man/L_OR.Rd
 7d39bbd12ade83dd6b9db384db5f5de8 *man/L_RR.Rd
 cffce60abba6c46e9e7148b1da01d5e8 *man/L_corr.Rd
 28e3830c1013a3d1b383d68191525fa7 *man/L_efficacy.Rd

NAMESPACE

@@ -16,6 +16,7 @@ export(L_t_test_sample_size)
 export(L_ttest)
 importFrom(graphics,curve)
 importFrom(graphics,lines)
+importFrom(graphics,plot)
 importFrom(graphics,segments)
 importFrom(stats,aggregate)
 importFrom(stats,anova)

R/L_1way_ANOVA.R

@@ -57,6 +57,7 @@
 #' @export
 #'
 #' @importFrom stats anova
+#' @importFrom graphics plot
 #' @importFrom stats lm
 #' @importFrom stats contr.poly
 #' @importFrom stats aggregate

R/L_1way_RM_ANOVA.R

@@ -48,6 +48,7 @@
 #' @export
 #'
 #' @importFrom stats anova
+#' @importFrom graphics plot
 #' @importFrom stats lm
 #' @importFrom stats aggregate
 #' @importFrom stats pf

R/L_1way_cat.R

@@ -1,7 +1,8 @@
 #' Likelihood Support for One-way Categorical Data
 #'
 #' This function calculates the support for one-way categorical data (multinomial), also
-#' gives chi-squared statistics. If there are only 2 categories then binomial information
+#' gives chi-squared and likelihood ratio test (G) statistics. If there are only 2
+#' categories then binomial information
 #' is given too with likelihood interval, including the likelihood-based % confidence
 #' interval. Support for the variance being more different than expected (Edwards p 187,
 #' Cahusac p 158) is also calculated.
@@ -33,6 +34,10 @@
 #'
 #' $p.value - p value for chi-squared.
 #'
+#' $LR.test = the likelihood ratio test statistic.
+#'
+#' $lrt.p = the p value for the likelihood ratio test statistic
+#'
 #' Additional outputs for binomial:
 #'
 #' $prob.val - MLE probability from data.
@@ -103,12 +108,16 @@ L_1way_cat <- function(obs, exp.p=NULL, L.int=2, alpha=0.05, toler=0.0001, verb=
   suppressWarnings(mc <- chisq.test(obs,p=exp.p,rescale.p = TRUE)) # suppress warnings
   p.value <- mc$p.value
   chi.s  <- unname(mc$statistic)
+  lrt <- 2*Sup  # likelihood ratio statistic
+  LRt_p <- 1-pchisq(lrt,df)
   toogood <- df/2*(log(df/chi.s)) - (df - chi.s)/2
   if (len!=2) {
     if(verb) cat("\nSupport for difference from expected, corrected for ", df, " df = ",  round(Supc,3), sep= "",
       "\n Support for variance differing more than expected = ",
       round(toogood,3), "\n\n Chi-square(", df, ") = ", chi.s,
-      ",  p = ", round(p.value,5), ", N = ", n, "\n ")
+      ",  p = ", round(p.value,5),
+      "\n Likelihood ratio test G(", df, ") = ", round(lrt,3),
+    ", p = ", round(LRt_p,5), ", N = ", n, "\n ")
   return(invisible(list(S.val = Supc, uncorrected.sup = Sup, df = df,
               observed = obs, exp.p = exp.p, too.good = toogood,
               chi.sq = chi.s, p.value = p.value)))
@@ -139,14 +148,16 @@ L_1way_cat <- function(obs, exp.p=NULL, L.int=2, alpha=0.05, toler=0.0001, verb=
       "\n Support for variance differing more than expected = ", round(toogood,3),
       "\n\n S-", L.int," likelihood interval (red line) from ",
       c(round(xmin1L$minimum,5), " to ", round(xmin2L$minimum,5)),
-      "\n\nChi-square(", df, ") = ", round(chi.s,3), ",  p = ", p.value,  ", N = ", n,
+      "\n\nChi-square(", df, ") = ", round(chi.s,3), ",  p = ", p.value,
+      "\n Likelihood ratio test G(", df, ") = ", round(lrt,3),
+      ", p = ", round(LRt_p,5),", N = ", n,
       "\n Likelihood-based ", 100*(1-alpha), "% confidence interval from ",
       c(round(xmin1$minimum,5), " to ", round(xmin2$minimum,5)), "\n ")
   invisible(list(S.val = Sup, uncorrected.sup = Sup, df = df, prob.val = p,
                  succ.fail = c(a,r), exp.p = exp.p,
                  like.int = c(xmin1L$minimum, xmin2L$minimum), like.int.spec = L.int,
                  too.good = toogood,
-                 chi.sq = chi.s, p.value = p.value,
+                 chi.sq = chi.s, p.value = p.value, LR.test = lrt, lrt.p = LRt_p,
                  conf.int = c(xmin1$minimum, xmin2$minimum), alpha.spec = alpha,
                  err.acc = c(xmin1$objective, xmin2$objective)))
 }

R/L_2way_Factorial_ANOVA.R

@@ -10,8 +10,8 @@
 #' Corrected support is given where appropriate, using
 #' Akaike's correction (Hurvich & Tsai (1989)). No correction is necessary for the
 #' two contrasts support since they both involve 1 parameter. Unequal group sizes are
-#' accommodated. F, p and partial eta-squared values are given for the two
-#' factors and their interaction.
+#' accommodated, using type III sums of squares. F, p and partial eta-squared values
+#' are given for the two factors and their interaction.
 #'
 #' @usage L_2way_Factorial_ANOVA(data, factor1, factor2, contrast1=NULL, contrast2=NULL, verb=TRUE)
 #'
@@ -118,28 +118,30 @@ L_2way_Factorial_ANOVA <- function(data, factor1, factor2, contrast1=NULL, contr
   F1 <- factor1
   F2 <- factor2
 
-  m1=anova(lm(dat~F1*F2))
+  options(contrasts = c("contr.sum","contr.poly")) # for type III SS
+  model <- lm(dat~F1*F2)
+  m1 <- drop1(model, .~., test="F")
 
-  tss <- sum(m1$`Sum Sq`)
-  eta_sq_1 <- m1$`Sum Sq`[1]/(m1$`Sum Sq`[1] + m1$`Sum Sq`[4]) # partial eta-squared
-  eta_sq_2 <- m1$`Sum Sq`[2]/(m1$`Sum Sq`[2] + m1$`Sum Sq`[4])
-  eta_sq_12 <- m1$`Sum Sq`[3]/(m1$`Sum Sq`[3] + m1$`Sum Sq`[4])
-  N <- (sum(m1$Df)+1)
-  ka <- m1$Df[1]+1
-  kb <- m1$Df[2]+1
+  tss <- sum(m1$`Sum of Sq`[2:4],m1$RSS[1])
+  eta_sq_1 <- m1$`Sum of Sq`[2]/(m1$`Sum of Sq`[2] + m1$RSS[1]) # partial eta-squared
+  eta_sq_2 <- m1$`Sum of Sq`[3]/(m1$`Sum of Sq`[3] + m1$RSS[1])
+  eta_sq_12 <- m1$`Sum of Sq`[4]/(m1$`Sum of Sq`[4] + m1$RSS[1])
+  N <- length(dat)
+  ka <- m1$Df[2]+1
+  kb <- m1$Df[3]+1
 
-  S_12 <- -0.5 * N * (log(m1$`Sum Sq`[4]) - log(tss))
+  S_12 <- -0.5 * N * (log(m1$RSS[1]) - log(tss))
   k2 <- 2       # parameters for variance and grand mean
-  k1 <- sum(m1$Df[1:3]) + k2
+  k1 <- sum(m1$Df[2:4]) + k2
 
   # Akaike's correction
   Ac <- function(k1,k2,N) { k2 * N/(N - k2 - 1) - k1 * (N/(N - k1 - 1)) }
 
   S_12c <- S_12 + Ac(k1,k2,N)
 
   # full model versus main effects only
-  S_FM <- -0.5 * N * (log(m1$`Sum Sq`[4]) - log(sum(m1$`Sum Sq`[3:4])))
-  k2 <- 2 + sum(m1$Df[1:2])
+  S_FM <- -0.5 * N * (log(m1$RSS[1]) - log(sum(m1$`Sum of Sq`[4],m1$RSS[1])))
+  k2 <- 2 + sum(m1$Df[2:3])
   S_FMc <- S_FM + Ac(k1,k2,N)
 
   datf <- data.frame(dat, F1, F2)
@@ -158,9 +160,9 @@ L_2way_Factorial_ANOVA <- function(data, factor1, factor2, contrast1=NULL, contr
     SS_cont1 <- sum(contrast1*means)^2/(sum(contrast1^2/(numbers)))
     SS_c1_unex <- tss - SS_cont1
 
-    S_c1M <- -0.5 * N * (log(SS_c1_unex) - log(sum(m1$`Sum Sq`[3:4])))
+    S_c1M <- -0.5 * N * (log(SS_c1_unex) - log(sum(m1$`Sum of Sq`[4],m1$RSS[1])))
     k1 <- 2 + 1    #2 parameters (variance and grand mean) + one for contrast
-    k2 <- 2 + sum(m1$Df[1:2])
+    k2 <- 2 + sum(m1$Df[2:3])
     S_c1Mc <- S_c1M + Ac(k1,k2,N)
 
     if (!is.null(contrast2)) {
@@ -173,28 +175,30 @@ L_2way_Factorial_ANOVA <- function(data, factor1, factor2, contrast1=NULL, contr
 
   interaction.plot(F1, F2, dat, ylab="Means")
 
-  Fval <- m1$`F value`
-  Pval <- m1$`Pr(>F)`
-  dfv <- m1$Df
+  Fval <- m1$`F value`[2:4]
+  Pval <- m1$`Pr(>F)`[2:4]
+  dfv <- m1$Df[2:4]
+  dfres <- (N-sum(dfv)-1)
 
-  Fval_c1 <- unname(SS_cont1/m1$`Mean Sq`[4])
-  Pval_c1 <- pf(Fval_c1, 1, m1$Df[4], lower.tail = FALSE)
+  res_msq <- m1$RSS[1]/dfres
+  Fval_c1 <- unname(SS_cont1/res_msq)
+  Pval_c1 <- pf(Fval_c1, 1, dfres, lower.tail = FALSE)
 
   if(verb) cat("\nSupport for full model (including interaction) versus null = ", round(S_12c,3), sep= "",
       "\n Support for full model versus main effects = ", round(S_FMc,3),
       "\n Support for contrast 1 versus main effects = ", if (!is.null(contrast1)) round(S_c1Mc,3),
       "\n Support for contrast 1 versus contrast 2 = ", if (!is.null(contrast2)) round(S_c1c2,3),
-      "\n\nFirst factor main effect F(", dfv[1],",", dfv[4],") = ", round(Fval[1],3),
+      "\n\nFirst factor main effect F(", dfv[1],",", dfres,") = ", round(Fval[1],3),
       ", p = ", Pval[1], ", partial eta-squared = ", round(eta_sq_1,3),
-      "\n Second factor main effect F(", dfv[2],",", dfv[4],") = ", round(Fval[2],3),
+      "\n Second factor main effect F(", dfv[2],",", dfres,") = ", round(Fval[2],3),
       ", p = ", Pval[2], ", partial eta-squared = ", round(eta_sq_2,3),
-      "\n Interaction F(", dfv[3],",", dfv[4],") = ", round(Fval[3],3),
+      "\n Interaction F(", dfv[3],",", dfres,") = ", round(Fval[3],3),
       ", p = ", Pval[3],  ", partial eta-squared = ", round(eta_sq_12,3),
-      "\n Contrast 1 F(1,",dfv[4],") = ", if (!is.null(contrast1)) round(Fval_c1,3),
+      "\n Contrast 1 F(1,",dfres,") = ", if (!is.null(contrast1)) round(Fval_c1,3),
       ", p = ", Pval_c1, "\n ")
 
   invisible(list(S.12c = S_12c, S.12 = S_12, S.FMc = S_FMc, S.FM = S_FM,
-                 S.c1.Mc = S_c1Mc, S.c1.M = S_c1M, S.c1.c2 = S_c1c2, Means = mean_out, df = m1$Df,
+                 S.c1.Mc = S_c1Mc, S.c1.M = S_c1M, S.c1.c2 = S_c1c2, Means = mean_out, df = c(m1$Df[2:4],dfres),
                  F.f1 = Fval[1], Pval.f1 = Pval[1], eta.sq.1 = eta_sq_1, F.f2 = Fval[2],
                  Pval.f2 = Pval[2], eta.sq.2 = eta_sq_2,
                  F.int = Fval[3], Pval.int = Pval[3], eta.sq.12 = eta_sq_12,

R/L_2way_cat.R

@@ -6,7 +6,7 @@
 #' is calculated. The support
 #' for trend across the columns is given (assuming the levels for columns are ordered),
 #' and conventional p value for trend.
-#' Finally, the usual chi-squared statistic and p value are given.
+#' Finally, Chi-squared and likelihood ratio test (G) statistics are given.
 #'
 #' @usage L_2way_cat(table, verb=TRUE)
 #'
@@ -44,9 +44,13 @@
 #'
 #' $residuals - the Pearson residuals.
 #'
-#' $chi.sq = the chi-squared statistic.
+#' $LR.test = the likelihood ratio test statistic.
 #'
-#' $p.value - the p value associated with the chi-squared statistic.
+#' $lrt.p - the p value for likelihood ratio test.
+#'
+#' $chi.sq - chi-squared value.
+#'
+#' $p.value - p value for chi-squared.
 #'
 #' $trend.p - p value for trend (from chi-squared dist.).
 #'
@@ -83,7 +87,7 @@ L_2way_cat <- function(table, verb=TRUE) {
 
 # calculating the interaction
   S2way <- 0
-  suppressWarnings(lt <- chisq.test(table)) # ignore warning message
+  suppressWarnings(lt <- chisq.test(table, correct=FALSE)) # ignore warning message
 # lt$observed * log(lt$observed/lt$expected) # individual S terms
   S2way <- sum( lt$observed * log(lt$observed/lt$expected) )
   df <- unname(lt$parameter)
@@ -108,28 +112,38 @@ L_2way_cat <- function(table, verb=TRUE) {
   chi.s <- unname(lt$statistic)
   toogood <- df/2*(log(df/chi.s)) - (df - chi.s)/2
 
+# likelihood ratio test
+  lrt <- 2*S2way
+  LRt_p <- 1-pchisq(lrt,df)
+
 # evidence for trend
-  table.pos <- table[1:length(col_sum)]
-  trX <- prop.trend.test(table[1,], col_sum)
-  tr <- unname(trX$statistic)/2      # S for trend
+  trX <- NULL
+  tr <- NULL
+  trX$p.value <- NULL
+  if (length(col_sum)>=3) {
+    table.pos <- table[1:length(col_sum)]
+    trX <- prop.trend.test(table[1,], col_sum)
+    tr <- unname(trX$statistic)/2      # S for trend
+  }
 
   if(verb) cat("\nSupport for interaction corrected for ", df, " df = ", round(S2wayc,3), sep= "",
     "\n Support for rows main effect corrected for ",
     length(row_sum)-1, " df = ", round(RowMain_c,3),
     "\n Support for columns main effect corrected for ", length(col_sum)-1, " df = ",
     round(ColMain_c,3), "\n Total support for whole table = ", round(Tot_S,3),
-    "\n Support for trend across columns = ", round(tr,3),
+    "\n Support for trend across columns = ", if (length(col_sum)>=3) round(tr,3),
     "\n Support for variance differing more than expected = ", round(toogood,3),
-    "\n\n P value for trend from chi-squared = ", trX$p.value,
-    "\n Chi-square(", df, ") = ", round(chi.s,3),
-    ",  p = ", round(lt$p.value,5), ", N = ", grandtot, "\n ")
+    "\n\n Chi-squared(", df, ") = ", round(chi.s,3), ",  p = ", round(lt$p.value,5),
+    "\n Likelihood ratio test G(", df, ") = ", round(lrt,3),
+    ", p = ", round(LRt_p,5), ", N = ", grandtot,
+    "\n\n Trend p value from chi-squared = ", if (length(col_sum)>=3) trX$p.value, "\n ")
 
   invisible(list(S.int = S2wayc, df = df, S.int.unc = S2way,
                S.Main.rows = RowMain_c, S.Main.cols = ColMain_c,
                S.Mr.uncorr = RowMain, S.Mc.uncorr = ColMain,
                df.rows = length(row_sum)-1, df.cols = length(col_sum)-1,
                S.total = Tot_S, S.trend = tr, too.good = toogood,
                observed = lt$observed, expected = lt$expected,
-               residuals = lt$residuals,
-                      chi.sq = chi.s, p.value = lt$p.value, trend.p = trX$p.value))
+               residuals = lt$residuals, chi.sq = lt$statistic, p.value = lt$p.value,
+               LR.test = lrt, lrt.p = LRt_p, trend.p = trX$p.value))
 }

R/L_OR.R

@@ -5,7 +5,7 @@
 #' and null (which is assumed to be 1, but can also be specified) values. A likelihood function
 #' is plotted for the obtained OR with a specified likelihood interval, and expected OR,
 #' if specified.
-#' Chi-squared statistics are also provided and a likelihood-based % confidence interval.
+#' Chi-squared and likelihood ratio test (G) statistics are also provided and a likelihood-based % confidence interval.
 #' It uses the optimize function to locate desired limits for both intervals and other
 #' support calculations.
 #'
@@ -45,6 +45,10 @@
 #'
 #' $p.value - p value.
 #'
+#' $LR.test = the likelihood ratio test statistic.
+#'
+#' $lrt.p = the p value for the likelihood ratio test statistic
+#'
 #' $residuals - the Pearson residuals.
 #'
 #' $alpha - specified significance level.
@@ -58,6 +62,7 @@
 #' @export
 #'
 #' @importFrom graphics segments
+#' @importFrom graphics plot
 #' @importFrom graphics lines
 #' @importFrom stats optimize
 #' @importFrom stats chisq.test
@@ -155,6 +160,8 @@ L_OR <- function(table, null=1, exp.OR=NULL, L.int=2, alpha=0.05, cc=FALSE, tole
   nullh <- exp(-sum(a*log(a/xa), b*log(b/(c1tot-xa)), c*log(c/(r1tot-xa)), d*log(d/(r2tot-c1tot+xa))))
 
   S2way <- log(1) - log(nullh) # check that this should be the same as S for observed OR
+  lrt <- 2*S2way  # likelihood ratio statistic
+  LRt_p <- 1-pchisq(lrt,1)
 
 # do the plot with lines
   plot(xs,ys,xlim=c(lolim,hilim),type="l", lwd = 1, xlab = "Odds Ratio", ylab = "Likelihood")
@@ -181,7 +188,9 @@ L_OR <- function(table, null=1, exp.OR=NULL, L.int=2, alpha=0.05, cc=FALSE, tole
       if (!is.null(exp.OR)) round(SexOR_null,3),
       sep= "", "\n   S-", L.int," likelihood interval (red line) is from ",
       c(round(begL,5), " to ", round(endL,5)),
-      "\n\nChi-square(1) = ", round(lt$statistic,3), ",  p = ", round(lt$p.value,4), ", N = ", grandtot,
+      "\n\nChi-square(1) = ", round(lt$statistic,3), ",  p = ", round(lt$p.value,4),
+      "\n Likelihood ratio test G(1) = ", round(lrt,3),
+      ", p = ", round(LRt_p,5),", N = ", grandtot,
       "\n   Likelihood-based ", 100*(1-alpha), "% confidence interval from ",
       c(round(beg,5), " to ", round(end,5)), "\n ")
 
@@ -190,6 +199,7 @@ L_OR <- function(table, null=1, exp.OR=NULL, L.int=2, alpha=0.05, cc=FALSE, tole
                S.exp.ORvsNull = SexOR_null, L.int = c(begL, endL), S_int = L.int,
                observed = lt$observed, expected = lt$expected,
                chi.sq = lt$statistic, corrected = cc, p.value = lt$p.value,
+               LR.test = lrt, lrt.p = LRt_p,
                residuals = lt$residuals, alpha = alpha, conf.int = c(beg, end),
                all.err.acc = c(xmin1$objective, xmin2$objective,
                                xmin1L$objective, xmin2L$objective,

R/L_RR.R

@@ -58,6 +58,7 @@
 #' @export
 #'
 #' @importFrom graphics segments
+#' @importFrom graphics plot
 #' @importFrom graphics lines
 #' @importFrom stats optimize
 #' @importFrom stats chisq.test