BiostatR Blog

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. In conclusion, RnetCDF does not scale data using scale_factor and add_offset automatically. Take care. 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 # take care, it does not scale automatically the data point_SST.RNetCDF_wrong <- var.get.nc(ncfile=dat.RNetCDF, variable="sst",                                  start=c(indice ,1), count=c(1, 1, 1)) # correct method point_SST.RNetCDF_corr <

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