BiostatR Blog

 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...

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 ...

Let use the same nc file obtained from NOAA ( http://www.esrl.noaa.gov/psd/data/gridded/data.noaa.oisst.v2.highres.html ) be used. Here I show the 3 ways to get annual data depending on the package used. utcal.nc from package RnetCDF is very useful. dest <- "sst.day.mean.2014.v2.nc" For all the package, I will use the same coordinates: lat=35, lon=28 # Note that the function ind_long_lat can be used data from any of these packages library(RNetCDF) dat.RNetCDF <- open.nc(dest) lat=35 lon=28 # it will be used for all data library(phenology) indice <- ind_long_lat(ncdf=dat.RNetCDF, long=lon, lat=lat) Packages RnetCDF point_SST.RNetCDF_corr <- var.get.nc(ncfile=dat.RNetCDF, variable="sst",                                       start=c(indice ,1), count=c(1, 1, 1)) # Get the time label.time <- dim.inq.nc(dat.RNetCDF, 0)$name date.char <- utcal.nc(att.get.nc(dat.RNet...

Just add \\ at the position of return \documentclass[a4paper]{article} \usepackage[utf8]{inputenc} \begin{document} \tableofcontents \section{Essai\\Pour un titre} <<>>= degree <- "o" print(degree) @ \end{document}

If you need fit a probability using optim(), it is more secure to fit an inverse logit value and use a logit transformation after to ensure that the values stay always between 0 and 1: L <- c(1.2) p <- 1/(1+exp(L)) And then the probabilities of the two events are p and 1-p. However it is more difficult in the case of multinomial distribution because the sum of the probabilities must be 1. Here is the correct way to do. For N events, you need N-1 parameters: L <- c(10, -1, 3, 7, 0) p <- 1/(1+exp(L)) p [1] 4.539787e-05 7.310586e-01 4.742587e-02 9.110512e-04 5.000000e-01 The trick is to transform these probabilities into conditional probabilities: p.multinomial <- function(p) c(p, 1)*c(1, sapply(seq_along(p), function(i) prod(1-p[1:i]))) p.multinomial(p) [1] 4.539787e-05 7.310254e-01 1.275420e-02 2.333885e-04 1.279708e-01 1.279708e-01 sum(c(p, 1)*c(1, sapply(seq_along(p), function(i) prod(1-p[1:i])))) [1] 1 To do the reverse: inv.p.multinomial <- f...

save.image(file= " xxx.Rdata") save the content of memory to be reused, but it does not work within a function. There, you can use: save(list = ls(all.names = TRUE), file = "total.Rdata", envir = environment()) For example: essai <- function(x) {   blabla <- c("ldklklkdlkd")   save(list = ls(all.names = TRUE), file = "total.Rdata", envir = environment()) } load("total.Rdata")

# generate data jo <- 1:800 trans_trigo <- 10+sin((jo/365.25)*2*pi)*1+ cos((jo/365.25)*2*pi)*5 obs <- rnorm(n=800, mean=trans_trigo, sd=2) # plot the data plot(jo, obs, bty="n", xlim=c(1,800), ylim=c(0,20)) par(new=TRUE) plot(jo, trans_trigo, bty="n", xlim=c(1,800), ylim=c(0,20), xlab="", ylab="",      axes=FALSE, col="violet", type="l", lwd=2) # group is per month observation mydat <- data.frame(obs=obs, jour=jo, jour2=jo*jo, group=round(jo/27)) # glm with ordinal day as linear factor mod <- glm(obs ~ jour, data=mydat, family=gaussian) summary(mod) newdata <- data.frame(jour=1:800) preddf2 <- predict(mod, type="response", newdata=newdata, se.fit=FALSE) par(new=TRUE) plot(jo, preddf2, bty="n", xlim=c(1,800),  ylim=c(0,20), xlab="", ylab="", axes=FALSE, col="red", type="l") # glm with ordinal day as second degree polynom mod3...