Rでプログラムでけた！！

　なかなか難産でしたが、Ｒで組んでたプログラムどうにか出来ました。

　機能としては、株の大きさからべき乗式使って胸高直径を推定。うんで、その胸高直径からさらに測っていない樹高を今度は拡張相対成長式使って推定するという合わせ技。しかもどっちも非線形回帰。非線形なのでPerlでもできるかもしれんが、多分めんどくさそう。といって、プロットごとに回帰させるのもちょっと。というわけで、Ｒで組んでみようということになったんですな～

　使うデータセットはDBH（間伐木の分は当然欠損）、根株の直径、樹高（欠損値あり）、それと各プロットのＩＤ。

　Ｒで組むときの大きなメリットは再利用のしやすさでしょうね。今回のやり方も全てデータフレームで管理するようにしてるので、例えば、ヘッダーのＩＤをそろえると（今回の場合はname DBH GBH,Hとか）他のデータでも応用が簡単にできます。そこは、基本的に１次元のデータ構造しか持たないPerlからすると、大きなメリットです。ま、その分データ構造で挙動が違ってくるので、そこの扱いに慣れるのは大変ですけど・・・

　それと、解析結果の二次利用が簡単ということでしょうか。今回の場合は非線形回帰で得られたパラメータを使って推定に使うというまさしく非線形回帰の結果を二次的に利用しているのですが、パラメータを取り出すのもcoef()で簡単に実現してしまいます。

　もしかすると、PerlでもCPANとかで探すと見つかるのかもしれませんが（ちなみに僕が探したときは見つかりませんでした）、マニアックな解析してそその結果を他の処理に渡す・・・なんてやつはＲプログラムの独壇場でしょう。

　ということで、書いたソースを追いときます。ちなみにヘッダーは以下のような感じです。

　　　　name,Subplot,色,No,sp,GBH,G 0.3,H,under.c,備考

　汚いので晒せるもんでもないですが、つっこみ大歓迎ということで・・・・
　

setwd(“d:”)
library(nls)
dat<-read.csv("testdata.csv",header=T) ##毎木データをプロットごとのデータフレームに分ける ############################################################# #name列はプロットＩＤ列。summaryの結果をname.matrixにぶっ込む name.matrix <- summary (dat$name) name.flame <- data.frame(name.matrix) ##プロットＩＤのそれぞれをベクトルにしてname.datにぶっ込む name.dat <- rownames(name.flame) ##各プロットの毎木データをそれぞれのデータフレームに分ける x <- 0 for( i in 1:length(name.dat)){ x <- x+1 each.name <- name.dat[x] each.dataframe <- subset(dat,dat$name == name.dat[x] ) ##assign()でプロットＩＤをオブジェクト名に当てる。 assign(each.name,each.dataframe) } ##一番外側のループ。各プロットデータに充てたデータフレームをループ ##今回は各プロットでべき乗式、拡張相対成長式を当てはめ、欠損値を補完する。 ############################ y <- 0 for( k in 1:length(name.dat)){ y <- y + 1 ##間伐している木の胸高を推定し、欠損値を補完 ##################################################### res.nls <- nls(GBH~A*G.0.3^B,data=get(name.dat[y]),start=c(A=1,B=1)) ##回帰係数をとる res.coef<-coef(res.nls) ##ラベルをとる res.coef <- as.vector(res.coef) par(ask=TRUE) plot(get(name.dat[y])$GBH~get(name.dat[y])$G.0.3) line <- function(x) {res.coef[1]*x^res.coef[2]} x <- seq(0,,length=100) lines(x,line(x),col="red") x <- 0 mod.GBH.vec<- NULL for( mod.GBH in 1:length(get(name.dat[y])$GBH)){ x <- x+1 if (is.na(get(name.dat[y])$GBH[x])==TRUE){ mod.GBH.vec <- append(mod.GBH.vec,res.coef[1]*get(name.dat[y])$G.0.3[x]^res.coef[2]) } else{ mod.GBH.vec <- append(mod.GBH.vec,get(name.dat[y])$GBH[x]) } } temp.data <- cbind(get(name.dat[y]),mod.GBH.vec) temp.data$DBH <- temp.data$mod.GBH.vec/pi assign(name.dat[y],temp.data) ##樹高未測定木の樹高を推定し、欠損値を補完 ##################################################### sub.forH <- subset(get(name.dat[y]),is.na(H)==FALSE) edit(sub.forH) inv.H <- 1/sub.forH$H resH.nls <- nls(inv.H ~ 1/(A * sub.forH$DBH^B)+ 1/C ,data = sub.forH,start=c(A=1,B=1,C=25)) resH.nls ##回帰係数をとる res.coef<-coef(resH.nls) ##ラベルをとる res.coef <- as.vector(res.coef) res.coef ##散布図と回帰線を描画 par(ask=TRUE) plot(sub.forH$H~sub.forH$DBH,xlim=c(0,50),ylim=c(0,30)) x <- seq(0,50,length=100) line <- function(x) {1/(1/(res.coef[1]*x^res.coef[2])+1/res.coef[3])} lines(x,line(x),col="red") x <- 0 mod.H <- NULL for( j in 1:length(get(name.dat[y])$name)){ x <- x+1 if (is.na(get(name.dat[y])$H[x])==TRUE){ mod.H <- append(mod.H,1/(1/(res.coef[1]*get(name.dat[y])$DBH[x]^res.coef[2])+1/res.coef[3])) } else{ mod.H <- append(mod.H,get(name.dat[y])$H[x]) } } temp.data <- cbind(get(name.dat[y]),mod.H) assign(name.dat[y],temp.data) edit(get(name.dat[y])) }