Bivariate confidence interval as a bubble (irregular ellipse)

This function will be useful if you want show the confidence interval of bivariate plots as a bubble (an irregular ellipse).
Note that this function has been introduced and enhanced in my package HelpersMG available in CRAN

ellipse <- function(center.x = 0, center.y = 0,
                        radius.x = 1, radius.y = 1,
                        radius.x.lower=NULL, radius.x.upper=NULL,
                        radius.y.lower=NULL, radius.y.upper=NULL, col=rgb(0, 0, 0, 0.1)) {

  if (is.null(radius.x.lower))  radius.x.lower <- radius.x
  if (is.null(radius.y.lower))  radius.y.lower <- radius.y
  if (is.null(radius.x.upper))  radius.x.upper <- radius.x
  if (is.null(radius.y.upper))  radius.y.upper <- radius.y

  for (k in seq_along(center.x)) {

  theta <- seq(0, pi / 2, length=100)
  x <- center.x[k] + radius.x.upper[k] * cos(theta)
  y <- center.y[k] + radius.y.upper[k] * sin(theta)

  theta <- seq(pi / 2, pi, length=100)
  x <- c(x, center.x[k] + radius.x.lower[k] * cos(theta))
  y <- c(y, center.y[k] + radius.y.upper[k] * sin(theta))

  theta <- seq(pi, 3/2*pi, length=100)
  x <- c(x, center.x[k] + radius.x.lower[k] * cos(theta))
  y <- c(y, center.y[k] + radius.y.lower[k] * sin(theta))

  theta <- seq(3/2*pi, 2 * pi, length=100)
  x <- c(x, center.x[k] + radius.x.upper[k] * cos(theta))
  y <- c(y, center.y[k] + radius.y.lower[k] * sin(theta))


  polygon(x, y, border=NA, col=col)
  }
}




Here are some examples

ellipse(center.x = c(0.2, 0.3, 0.25), center.y = c(0.7, 0.6, 0.55), radius.x = c(0.1, 0.1, 0.1), radius.y = c(0.15, 0.2, 0.5))

ellipse(center.x = 0.5, center.y = 0.5, radius.x.lower = 0.1, radius.x.upper = 0.3, radius.y = 0.2)
ellipse(center.x = 0.6, center.y = 0.3, radius.x.lower = 0.1, radius.x.upper = 0.3, radius.y.lower = 0.2, radius.y.upper = 0.4)


males <- c(10, 25, 3, 4)
N <- c(12, 52, 17, 10)


males2 <- c(12, 20, 3, 6)
N2 <- c(15, 50, 20, 12)

bc <- HelpersMG::.BinomialConfidence(x=males, n=N)
bc2 <- HelpersMG::.BinomialConfidence(x=males2, n=N2)

plot(0:1, 0:1, xlim=c(0, 1), ylim=c(0,1), lty=2, type="l", las=1, bty="n")

ellipse(center.x = bc[, "PointEst"], center.y = bc2[, "PointEst"],
            radius.x.lower = bc[, "PointEst"]-bc[, "Lower"],
            radius.x.upper = bc[, "Upper"]-bc[, "PointEst"],
            radius.y.lower = bc2[, "PointEst"]-bc2[, "Lower"],
            radius.y.upper = bc2[, "Upper"]-bc2[, "PointEst"])

Commentaires

Enregistrer un commentaire

Posts les plus consultés de ce blog

Standard error from Hessian Matrix... what can be done when problem occurs

An Hessian matrix is a square matrix of partial second order derivative. From the square-root of the inverse of the diagonal, it is possible to estimate the standard error of parameters. Let take an example: val=rnorm(100, mean=20, sd=5) # Return -ln L of values in val in Gaussian distribution with mean and sd in par fitnorm<-function(par, val) {   -sum(dnorm(val, par["mean"], par["sd"], log = TRUE)) } # Initial values for search p<-c(mean=20, sd=5) # fit the model result <- optim(par=p, fn=fitnorm, val=val, method="BFGS", hessian=TRUE) # Inverse the hessian matrix to get SE for each parameters mathessian <- result$hessian inversemathessian <- solve(mathessian) res <- sqrt(diag(inversemathessian)) # results data.frame(Mean=result$par, SE=res)   library(MASS) (r<-fitdistr(val, "normal")) It works well. However, the inverse of the matrix cannot be calculated if the second order derivative for some paramete
Lire la suite

stepAIC from package MASS with AICc

Due to a bug in dropterm() and addterm() in package MASS, it is impossible to use AICc with stepAIC() . Here is a solution. Enter these commands and now you can use: stepAIC(g0, criteria="AICc") or stepAIC(g0, criteria="BIC") For example: datax <- data.frame(y=rnorm(100), x1=rnorm(100),                     x2=rnorm(100), x3=rnorm(100), x4=rnorm(100)) g0 <- glm(y ~ x1+x2+x3+x4, data=datax) gx_AICc <- stepAIC(g0, criteria="AICc") gx_AIC <- stepAIC(g0, criteria="AIC") gx_AIC <- stepAIC(g0, criteria="BIC")
Lire la suite

Install treemix in ubuntu 20.04

Informations are available here  https://bitbucket.org/nygcresearch/treemix/downloads/ In terminal, enter for  Ubuntu 20.04: wget https://bitbucket.org/nygcresearch/treemix/downloads/treemix-1.13.tar.gz tar -xvf treemix-1.13.tar.gz sudo apt-get install libgsl-dev sudo apt install libboost-dev sudo apt install libboost-all-dev cd treemix-1.13/ ./configure make sudo make install And it works: $ treemix TreeMix v. 1.13 $Revision: 231 $ Options: -h display this help -i [file name] input file -o [stem] output stem (will be [stem].treeout.gz, [stem].cov.gz, [stem].modelcov.gz) -k [int] number of SNPs per block for estimation of covariance matrix (1) -global Do a round of global rearrangements after adding all populations -tf [file name] Read the tree topology from a file, rather than estimating it -m [int] number of migration edges to add (0) -root [string] comma-delimited list of populations to set on one side of the root (for migration) -g [vertices file name] [edges file name] read the g
Lire la suite