R variable not found, but specifically defined

Question

I have written a function to run phylogenetic generalized least squares, and everything looks like it should work fine, but for some reason, a specific variable which is defined in the script (W) keeps coming up as undefined. I have stared at this code for hours and cannot figure out where the problem is.

Any ideas?

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)

returns "Error in eval(expr, envir, enclos) : object 'W' not found"

even though W is specifically defined!

Thanks!

Debugging code that one can't run on one's own computer is a special circle of hell that I generally aim to avoid. — joran
– joran, Commented May 16, 2012 at 23:42
You probably want cbind(dat, W) in the corMartins2() function. I think gls() is looking for W in the dat data.frame and not finding it. — jthetzel
– jthetzel, Commented May 16, 2012 at 23:45
sorry about that, I removed lowerB and upperB and replaced them with hard-coded values. Also, this should be a runable script. — Pascal
– Pascal, Commented May 16, 2012 at 23:48

mathematical.coffee · Accepted Answer · 2012-05-16 23:50:19Z

6

The error's occuring in the gls calls in myou and corMatrins2: you have to pass in W as a column in dat because gls is looking for it there (when you put weights = ~W as a formula like that it looks for dat$W and can't find it).

Just change data=dat to data=cbind(dat,W=W) in both functions.

answered May 16, 2012 at 23:50

mathematical.coffee

57.1k15 gold badges160 silver badges197 bronze badges

Sign up to request clarification or add additional context in comments.

1 Comment

Pascal Over a year ago

Thanks for the fix! This makes complete sense!

jthetzel · Accepted Answer · 2012-05-16 23:52:32Z

5

The example is not reproducible for me, as lowerB and upperB are not defined, however, perhaps the following will work for you, cbinding dat with 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)

answered May 16, 2012 at 23:52

jthetzel

3,6233 gold badges27 silver badges38 bronze badges

1 Comment

Pascal Over a year ago

yeah that was the problem...thanks for taking the time to help me!

Collectives™ on Stack Overflow

R variable not found, but specifically defined

2 Answers 2

1 Comment

1 Comment

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

1 Comment

1 Comment

Your Answer

Sign up or log in

Post as a guest

Related