Errors running R scripts

[r.a....@exeter.ac.uk]

Dear Joao

 

I am in need of some help with running two of the R scripts for popABC.

 

When running model_choice.r and get_modes.r I get the following errors:

 

 

> #Bar plot for model-choice (regression)

>      res.topol<-calmod(target,abc.rej[,2],abc.rej[,23:52],1,rej=F)

Error in calmod(target, abc.rej[, 2], abc.rej[, 23:52], 1, rej = F) :

  could not find function "predict.vglm"

>      print(res.topol$x2)

Error in print(res.topol$x2) : object 'res.topol' not found

>      write(res.topol$x2,"modelprob_reg.txt")

Error in cat(x, file = file, sep = c(rep.int(sep, ncolumns - 1), "\n"),  :

  object 'res.topol' not found

>      barplot(res.topol$x2,names=c("Model1","Model2"))

Error in barplot(res.topol$x2, names = c("Model1", "Model2")) :

  object 'res.topol' not found

>      abline(h=1.0/2,col="red")

>      print("model-choice done.")

 

 

> # Plot line for splitting time and save it in a .eps file

>      abc.reg.tev <- makepd4(target,abc.rej[,10],abc.rej[,16:33],tol=1,rej=F,transf="none",bb=c(mintev,maxtev))

>     write("splitting time (tev):", "estimates.txt", append=T)

>      write(loc1statsx(x=abc.reg.tev$x,wt=abc.reg.tev$wt, prob=credible, alpha=1.0, xlim=c(mintev,maxtev)), "estimates.txt", append=T)

Error in loc1statsx(x = abc.reg.tev$x, wt = abc.reg.tev$wt, prob = credible,  :

  could not find function "predict.locfit"

>      print("tev done.")

[1] "tev done."

 

The rejection, priors and target files all load. With get_modes, the rejection function works, the problem is with the regression code.

 

It would be great if you can suggest how I can fix these errors.

 

Many Thanks

Andy

[Joao Sollari Lopes]

I think you just need to install and load the locfit library.
Joao

[r.a....@exeter.ac.uk]

Joao

 

I have already installed the locfit library.

 

Best Wishes

 

Andy

[pconverse...@gmail.com]

Andy,

I was having a similar problem. I have solved it by slightly modifying an R script the loc1statsx function was calling upon. Perhaps Joao can double check to see if what I did is acceptable.

Open loc2plot_d.r in TextWrangler or another text editing software. Search for the string "predict.locfit" and it should jump to areas in the code where the function is. Here is an example:

{

sc1 <- sqrt(var(x))

sc2 <- sqrt(var(y))

if(missing(maxk)) maxk <- 100

if(missing(xlim)) fit <- locfit(~x+y,alpha=alpha,scale=c(sc1,sc2),

maxk=maxk,mint=100,cut=0.8,maxit=100)

else fit <- locfit(~x+y,alpha=alpha,scale=c(sc1,sc2),

xlim=xlim,maxk=maxk,mint=100,cut=0.8,maxit=100)

# d1 <- (x-px)^2+(y-py)^2

# best <- d1 == min(d1)

# lev <- mean(fitted(fit)[best])

lev <- predict.locfit(fit,list(px,py))

slev <- sort(fitted(fit))

indic <- slev <= lev

sum(indic)/length(x)

}

Simply replace this string with "lev <- predict(fit,list(px,py))". The modified code should look like this:

{

sc1 <- sqrt(var(x))

sc2 <- sqrt(var(y))

if(missing(maxk)) maxk <- 100

if(missing(xlim)) fit <- locfit(~x+y,alpha=alpha,scale=c(sc1,sc2),

maxk=maxk,mint=100,cut=0.8,maxit=100)

else fit <- locfit(~x+y,alpha=alpha,scale=c(sc1,sc2),

xlim=xlim,maxk=maxk,mint=100,cut=0.8,maxit=100)

# d1 <- (x-px)^2+(y-py)^2

# best <- d1 == min(d1)

# lev <- mean(fitted(fit)[best])

lev <- predict(fit,list(px,py))

slev <- sort(fitted(fit))

indic <- slev <= lev

sum(indic)/length(x)

}

You'll need to replace every instance of predict.locfit with the modification above.

Hope this helps,

Cheers,

Paul

[Joao Sollari Lopes]

Thanks for that Paul.
You are absolutely right. You should change "predict.locfit" for "predict".
When loc2plot_d.r was written by Mark, both uses were fine. But now the correct use is "predict".
Sorry for that, Andy.

Another suggestion is to use the package "abc", its syntax was partially based on Mark Beaumont's scripts so it should deal with popABC files very easily.
Some differences in the results are to be expected since there are some details that are different between the approaches, in particular:
a) density distribution are obtain with function "density" instead of using "locfit" package;
b) standardization of the summary statistics are performed by dividing by the MAD instead of by the standard deviation.

On the plus side, package "abc" allows to perform the regression-step using a more sophisticated neuronal network regression. The other advantage is that it is being continually mainted as an R package available in CRAN.

Best,
Joao