Periodic GLM

# 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 <- glm(obs ~ jour+jour2, data=mydat, family=gaussian)
summary(mod3)
newdata <- data.frame(jour=1:800, jour2=(1:800)^2)
preddf2 <- predict(mod3, 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="green", type="l")

# glm with ordinal day transformed as period factor
mod2 <- glm(obs ~ sin((jour/365.25)*2*pi)+ cos((jour/365.25)*2*pi), data=mydat, family=gaussian)
summary(mod2)
newdata <- data.frame(jour=1:800)
preddf2 <- predict(mod2, 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="blue", type="l")

# Same with glmm
library("MASS")
mod4 <- glmmPQL(obs ~ sin((jour/365.25)*2*3.141593)+ cos((jour/365.25)*2*3.141593), random= ~ 1 | group, data=mydat, family=gaussian)
summary(mod4)
newdata <- data.frame(jour=1:800, group=round(jo/27))
preddf2 <- predict(mod4, 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="yellow", type="l")

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