2

系統発生の一般化最小二乗法を実行する関数を作成しましたが、すべて正常に動作するように見えますが、何らかの理由で、スクリプト (W) で定義されている特定の変数が未定義として表示され続けます。私はこのコードを何時間も見つめてきましたが、どこに問題があるのか​​ わかりません.

何か案は？

myou <- function(alpha, datax, datay, tree){
    data.frame(datax[tree$tip.label,],datay[tree$tip.label,],row.names=tree$tip.label)->dat
    colnames(dat)<-c("Trait1","Trait2")
    W<-diag(vcv.phylo(tree)) # Weights
    fm <- gls(Trait1 ~ Trait2, data=dat, correlation = corMartins(alpha, tree, fixed = TRUE),weights = ~ W,method = "REML")
    return(as.numeric(fm$logLik))
}

corMartins2<-function(datax, datay, tree){
    data.frame(datax[tree$tip.label,],datay[tree$tip.label,],row.names=tree$tip.label)->dat
    colnames(dat)<-c("Trait1","Trait2")
    result <- optimize(f = myou, interval = c(0, 4), datax=datax,datay=datay, tree = tree, maximum = TRUE)
    W<-diag(vcv.phylo(tree)) # Weights
    fm <- gls(Trait1 ~ Trait2, data = dat, correlation = corMartins(result$maximum, tree, fixed =T),weights = ~ W,method = "REML")
    list(fm, result$maximum)}




#test


require(nlme)
require(phytools)
simtree<-rcoal(50)
as.data.frame(fastBM(simtree))->dat1
as.data.frame(fastBM(simtree))->dat2

corMartins2(dat1,dat2,tree=simtree)

「eval（expr、envir、enclos）のエラー：オブジェクト 'W'が見つかりません」を返します

W は具体的に定義されていますが!

ありがとう！

4

2 に答える 2

6

のgls呼び出しでエラーが発生しmyouている：そこで検索しているため、列としてcorMatrins2渡す必要があります（そのような数式として入力すると、検索して見つけることができません）。Wdatglsweights = ~Wdat$W

両方の機能でに変更data=datするだけです。data=cbind(dat,W=W)

于 2012-05-16T23:50:19.383 に答える
5

lowerB定義されていないため、この例は私には再現できませんupperBが、おそらく次のように機能しcbinding datますW。

myou <- function(alpha, datax, datay, tree){
    data.frame(datax[tree$tip.label,],datay[tree$tip.label,],row.names=tree$tip.label)->dat
    colnames(dat)<-c("Trait1","Trait2")
    W<-diag(vcv.phylo(tree)) # Weights
            ### cbind W to dat
            dat <- cbind(dat, W = W)
    fm <- gls(Trait1 ~ Trait2, data=dat, correlation = corMartins(alpha, tree, fixed = TRUE),weights = ~ W,method = "REML")
    return(as.numeric(fm$logLik))
}

corMartins2<-function(datax, datay, tree){
    data.frame(datax[tree$tip.label,],datay[tree$tip.label,],row.names=tree$tip.label)->dat
    colnames(dat)<-c("Trait1","Trait2")
    result <- optimize(f = myou, interval = c(lowerB, upperB), datax=datax,datay=datay, tree = tree, maximum = TRUE)
    W<-diag(vcv.phylo(tree)) # Weights
    ### cbind W to dat
    dat <- cbind(dat, W = W)
    fm <- gls(Trait1 ~ Trait2, data = dat, correlation = corMartins(result$maximum, tree, fixed =T),weights = ~ W,method = "REML")
    list(fm, result$maximum)}




#test
require(phytools)
simtree<-rcoal(50)
as.data.frame(fastBM(simtree))->dat1
as.data.frame(fastBM(simtree))->dat2

corMartins2(dat1,dat2,tree=simtree)
于 2012-05-16T23:52:32.353 に答える