Articles

Install tiff, ijtiff and jpeg packages in Ubuntu and MacOSX

In terminal, first install wget to download sources of the packages. Check the most recent sources in CRAN web site. brew install wget For tiff packages, in MacOSX, Intel computer  brew install libtiff sudo ln -s  /usr/local/include/tiff.h /Library/Frameworks/R.framework/Resources/include/ sudo ln -s /usr/local/include/tiffconf.h /Library/Frameworks/R.framework/Resources/include/ sudo ln -s /usr/local/include/tiffio.h /Library/Frameworks/R.framework/Resources/include/ sudo ln -s /usr/local/include/tiffvers.h /Library/Frameworks/R.framework/Resources/include/ sudo ln -s /usr/local/lib/libtiff* /Library/Frameworks/R.framework/Resources/lib/ cd $HOME wget https://cran.r-project.org/src/contrib/tiff_0.1-8.tar.gz R CMD INSTALL --configure-vars='INCLUDE_DIR=/usr/local/include LIB_DIR=/usr/local/lib' tiff_0.1-8.tar.gz rm tiff_0.1-8.tar.gz For ijtiff In R: install.packages(c("checkmate", "strex", "zeallot")) In terminal:  wget https://cran....

Setup new installation of R in MacosX

Adapted from  https://ryanhomer.github.io/posts/build-openmp-macos-catalina-complete Install R from  cran.r-project.org At the end of ~/.zshrc add:  export PATH=$PATH:/Library/Frameworks/R.framework/Versions/Current/Resources/bin/ Then you can use R within the terminal. Install macbrew via: /bin/bash -c "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/HEAD/install.sh)" Then in ~/.Rprofile: marcgirondot@MacBook-Pro-de-Marc ~ % cat .Rprofile # Default CRAN mirror options(repos=structure(c(CRAN="https://cloud.r-project.org/"))) # Path to run R with system Sys.setenv(PATH=paste(Sys.getenv("PATH"), "/Library/Frameworks/R.framework/Versions/Current/Resources/bin/", sep=":")) Install llvm brew install llvm libomp Install GCC and gettext brew install gcc gettext To configure R to build packages with the versions of clang and gcc you just installed, set up a Makevars file. This must be located at ~/.r/Makevars. Some parts do no w...

Generate a pseudo-Hessian matrix from MCMC series

The Hessian matrix at a point P (P is a vector) is the matrix of partial derivatives at this point. It is often used after non-linear optimization and the standard error of parameters can be obtained using the square-root of the diagonal of the inverse of the Hessian matrix. But when you are estimate SE like this, you lost the covariances between parameters. Then if you estimate random numbers using the SE of parameters, the original structure is lost. It is much better to estimate random numbers directly from the Hessian matrix. Then you keep the variance-covariance structure. When the model is fitted using MCMC, no Hessian matrix is directly available. A pseudo-Hessian matrix can be calculated using: hessian <- solve(cov( [matrix of values for different iterations] )) Let do an example with the fit of a normal distribution: First using maximum likelihood: val <- rnorm(30, 10, 2) library(MASS) ft <- fitdistr(val, densfun="normal") The Hessian matrix is not available...

Confidence interval to check for difference: Use 1.96 SE and not 1.96 SD !

Let create two variables of length 100, one with mean 10 (A) and one with mean 12 (B) both with SD=2. Of course the two variables overlap. A <- rnorm(100, 10, 2) B <- rnorm(100, 12, 2) library(HelpersMG) barplot_errbar(c(mean(A) ,mean(B)), errbar.y = c(1.96*sd(A), 1.96*sd(B)), las=1, ylim=c(0, 20), main="1.96 x SD") Now do a t test. It is highly significant: t.test(A, B) Welch Two Sample t-test data:  A and B t = -7.8344, df = 197.93, p-value = 2.82e-13 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval:  -2.737632 -1.636588 sample estimates: mean of x mean of y   9.741912 11.929021  Then if you want use the overall of confidence interval, use SE: barplot_errbar(c(mean(A) ,mean(B)), errbar.y = c(1.96*sd(A)/sqrt(length(A)),   1.96*sd(B)/sqrt(length(A))), las=1, ylim=c(0, 12), main="1.96 x SE")

What are the consequences of replacing missing data with median?

The conclusion is that it artificially reduced the variability of the correlation coefficient. It is bad practice. But it is much better than doing nothing ! cor.original <- NULL cor.na <- NULL cor.median <- NULL for (i in 1:10000) {   A <- rnorm(100, mean=100, sd=20)   B <- rnorm(100, mean=100, sd=20)   Bprime <- ifelse(sample(c(0,1), 100, replace = TRUE), B, NA)   Bter <- ifelse(is.na(Bprime), median(B, na.rm = TRUE), Bprime)   cor.original <- c(cor.original, cor(x=A, y=B, method = "spearman"))   cor.na <- c(cor.na, cor(x=A, y=Bprime, method = "spearman", use="complete.obs"))   cor.median <- c(cor.median, cor(x=A, y=Bter, method = "spearman", use="complete.obs")) } layout(1:3) hist(cor.original, xlim=c(-0.6, 0.6), breaks=seq(from=-0.6, to=0.6, by=0.05)) hist(cor.na, xlim=c(-0.6, 0.6), breaks=seq(from=-0.6, to=0.6, by=0.05)) hist(cor.median, xlim=c(-0.6, 0.6), breaks=seq(from=-0.6, to=0.6, by=0.05)) quanti...

More on linear regression... take care of initial point for robust regression

 It is better to supply a initial point for robust regression ! x <- c(0.428571428571429, 0.2, 0.3, 0, 0, 0.2, 0, 0, 0.1, 0, 0.1,    0.1, 0, 0, 0, 0.3, 0.2, 0.3, 0.2, 0.2, 0, 0, 0, 0.2, 0.222222222222222,    0.1, 0, 0.4, 0.3, 0.5, 0, 0.4, 0.5, 0.8, 0.3, 0.1, 0.2, 0, 0.1,    0, 0.1, 0.1, 0.4, 0, 0, 0, 0, 0, 0, 0.333333333333333, 0.444444444444444,    0.2, 0.222222222222222, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,    0, 0, 0, 0, 0.4, 0, 0, 0, 0, 0.111111111111111, 0, 0, 0, 0, 0,    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,    0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0.2, 0, 0, 0, 0, 0, 0, 0) y <- c(0.369999388816617, 0.152260672623962, 0.392636473518975, 0.107742543910461,         0.105802942749025, 0.147888875829182, 0.0180958542177892, 0.00376865991773073,         0.008...

Create a named vector with names stored in a variable

Often I want create a named vector and the names are stored in a variable. I must use a 2 steps function: names <- paste0("V", 1:100) A <- 1:100 names(A) <- names I search a way to do the same in a single step. Here is the solution: A <- structure(1:100, .Dim = 100L,                     .Dimnames = list(names)) But which one is the faster? names <- paste0("V", 1:100) library(microbenchmark) microbenchmark({     A <- 1:100     names(A) <- names},      {A <- structure(1:100, .Dim = 100L,                     .Dimnames = list(names))},      times = 1000L) Unit: microseconds                                                                 expr  min ...