Rで時系列データをcosカーブにフィッティングする

  • Rをインストールする。
    • ここ参考に。2015年現在たくさんネット上に情報があります。
    • RStudioが便利です。
  • 次にデータをタブ区切りのテキストファイルとしてこんな感じに用意します。エクセルで作ると楽ちんです。ファイルはお好きなところに置けばよいです。名前はとりあえず、test.txtとでもしておけば良いです。
  • 以下のソースコードをコンソールにコピー&ペースト
    • workingDirとinputfilenameのパスを書き換えましょう。
workingDir="/Users/hito/Dropbox/misc/r/" #ここは手元の環境にあわせて変えて下さい。
inputfilename="test.txt" #ここは手元の環境にあわせて変えて下さい。
setwd(workingDir)
fittingSin <- function(data,index,PDF){
dataNum=colnames(data)[index+1]
filename=paste(dataNum,".pdf",sep="")
print(filename)
x <- data$time
y <- data[,index+1]
x <- x[!is.na(y)]
y <- y[!is.na(y)]
xaxsmin <- ceiling(min(x)/12)
xaxsmax <- floor(max(x)/12)
xaxsvec <- xaxsmin:xaxsmax*12
result <- nls(y ~A*cos(2*pi/T*(x-a))+b, algorithm ="port", start=c(A=(max(y)-min(y))/2,a=2,T=24,b=mean(y)))
result1 <- try(nls(y ~A*cos(2*pi/T*(x-a))+b, algorithm ="port", start=c(A=(max(y)-min(y))/2,a=-6,T=24,b=mean(y))))
result2 <- try(nls(y ~A*cos(2*pi/T*(x-a))+b, algorithm ="port", start=c(A=(max(y)-min(y))/2,a=0,T=24,b=mean(y))))
result3 <- try(nls(y ~A*cos(2*pi/T*(x-a))+b, algorithm ="port", start=c(A=(max(y)-min(y))/2,a=6,T=24,b=mean(y))))
result4 <- try(nls(y ~A*cos(2*pi/T*(x-a))+b, algorithm ="port", start=c(A=(max(y)-min(y))/2,a=12,T=24,b=mean(y))))
result <- result1
if(deviance(result2)<deviance(result1)) {result <- result2}
if(deviance(result3)<deviance(result1)) {result <- result3}
if(deviance(result4)<deviance(result1)) {result <- result4} 
par = t(c(coef(result),deviance(result)))
rownames(par) = colnames(data[index+1]) 
parnames=c(names(coef(result)),"R^2")
predictdata <- predict(result)
if(PDF){pdf(filename)}
plot(x,y,pch=16,xlim=c(min(x),max(x)),ylim=c(0,0.6),yaxt="n",xaxt="n",xlab="Time[h]",ylab="Ratio of  P-KaiC",col="grey60")
axis(1,at=xaxsmin:xaxsmax*12,tck=1)
par(new=T)
plot(x,predictdata,col="red",ann=F,type="l",lwd=2,xaxt="n",xlim=c(min(x),max(x)),ylim=c(0,0.6))
sampleNo=sub("[^[:digit:]].","",colnames(data)[index+1])
mtext(paste("No.",sampleNo),cex=0.8)
textx <- (max(x)-min(x))*(-0.05) + min(x)
texty <- -(max(y)-min(y))*0.05 + 0.6
textpx <- (max(x)-min(x))*0.15 + min(x)
textpy <- -(max(y)-min(y))*0.05 + 0.6
options(digits=4)
fittext<- "y=A*cos(2*Pi/T*(t-a))+b\n"
textA <- paste("A=",round(par[1],digits=4),"\n")
texta <- paste("a=",round(par[2],digits=2),"\n")
textT <- paste("T=",round(par[3],digits=2),"\n")
textb <- paste("b=",round(par[4],digits=2),"\n")
text(textx,texty,fittext,pos=4)
text(textpx,textpy,paste(textA,texta,textT,textb),pos=1)
print("A,a,T,b")
print(c(round(par[1],digits=4),round(par[2],digits=2),round(par[3],digits=2),round(par[4],digits=2)))
if(PDF){dev.off()}
}

#get 1PDF
pdf("Fitting.pdf",family="Helvetica")
par(mfrow = c(5,4))
par(oma=c(0,0,0,0))
par(ann=FALSE)
par(mai=c(0.15,0,0.15,0))
par(mgp=c(1,0,0))
fitdata <- read.delim(inputfilename,header=T)
dataNum=length(colnames(fitdata))-1
for(i in 1:dataNum){
  PDF=FALSE
  fittingSin(fitdata,i,PDF)
}
dev.off()
  • PDFファイルができあがります。

Periodicな遺伝子発現のカーブを抽出してくれるRのライブラリ

ここ参考に


  編集 添付 リロード   新規 一覧 単語検索   最終更新のRSS