4

This is a very strange situation that I've come across. Basically, I'm trying to fit a cumulative distribution function to the G function of my data. After doing so, I want to plot the the model and the original data, and output this as PDF. I'll allow the code to explain (simply copy and paste):

library(spatstat)

data(swedishpines)

mydata <- swedishpines

mydata.Gest <- Gest(mydata)

Gvalues <- mydata.Gest$rs

count <- (which(Gvalues == 1))[1]

new_r <- seq(1/count, length(Gvalues)/count, by = 1/count)

GvsR_dataframe <- data.frame(G <- Gvalues, R <- new_r)

themodel <- suppressWarnings(nls(G ~ pnorm(R, mean, sd), data = GvsR_dataframe, start = list(mean=0.4, sd=0.2), trace = FALSE))

pdf(file = "ModelPlot.pdf")

plot(mydata.Gest, cbind(rs, theo) ~ new_r, lty = c(1, 2), col = c("black", "red"), xlim = c(0, max(new_r)), ylim = c(0,1), main = paste("Model-fitting for G Function \n Mean = ",as.numeric(coef(themodel)[1]),"\n Standard Deviation = ",as.numeric(coef(themodel)[2]), sep=''), ylab = "G(r)", xlab = "Distance Between Particles (r)", legend = NULL)
lines(new_r, predict(themodel), lty = 2, col = "blue")
legend("bottomright", c("CSR", "Swedish Pines", "Normal Probability \n Density Function"), lty = c(2, 4, 1, 2), col = c("red", "black", "blue"), bg = 'grey', border = 'black')

graphics.off()

The above code works perfectly.

Now for the strange part.

When I encapsulate all of the commands after mydata <- swedishpines as a function, and cause mydata to be an input to this function, it doesn't work any longer. The following code should perform just as the last segment of code did, but it does not.

library(spatstat)

data(swedishpines)

mydata <- swedishpines

ModelFit <- function(mydata) {

mydata.Gest <- Gest(mydata)

Gvalues <- mydata.Gest$rs

count <- (which(Gvalues == 1))[1]

new_r <- seq(1/count, length(Gvalues)/count, by = 1/count)

GvsR_dataframe <- data.frame(G <- Gvalues, R <- new_r)

themodel <- suppressWarnings(nls(G ~ pnorm(R, mean, sd), data = GvsR_dataframe, start = list(mean=0.4, sd=0.2), trace = FALSE))

pdf(file = "ModelPlot.pdf")

plot(mydata.Gest, cbind(rs, theo) ~ new_r, lty = c(1, 2), col = c("black", "red"), xlim = c(0, max(new_r)), ylim = c(0,1), main = paste("Model-fitting for G Function \n Mean = ",as.numeric(coef(themodel)[1]),"\n Standard Deviation = ",as.numeric(coef(themodel)[2]), sep=''), ylab = "G(r)", xlab = "Distance Between Particles (r)", legend = NULL)
lines(new_r, predict(themodel), lty = 2, col = "blue")
legend("bottomright", c("CSR", "Swedish Pines", "Normal Probability \n Density Function"), lty = c(2, 4, 1, 2), col = c("red", "black", "blue"), bg = 'grey', border = 'black')

graphics.off()

}

ModelFit(mydata)

The following error occurs:

Error in eval(expr, envir, enclos) : object 'new_r' not found

I'm VERY confused. I've been working on this for a long time, and just could not come up with a solution to this problem. The PDF is outputted, but it is corrupt, and will not open. I have no idea why new_r 'disappears', but in doing so, it causes all of the plotting operations to halt. Obviously new_r is local to the function ModelFit, but it almost seems as though it's local to certain areas in the function, as well.

Any help would be greatly appreciated.

MikeZ
  • 345
  • 1
  • 4
  • 15
  • I cannot reproduce. Are you certain the code in your second block is precisely the same code you are executing in R? – David Robinson Aug 02 '12 at 19:34
  • I am sure that it is. All I've changed in the second block is a function definition (which now contains the majority of the code from the first block), followed by a call to that function. – MikeZ Aug 02 '12 at 19:38
  • Interesting. I'd appreciate any help, as I'm trying to run a large program in batch mode on a linux server, and this small segment of code has caused me a lot of grief. – MikeZ Aug 02 '12 at 19:40
  • Disregard the above comments- I discovered that it does not throw the error only when I already had run the previous, working code in the same space. – David Robinson Aug 02 '12 at 19:40
  • Justin's answer below, putting in `mydata.Gest$r <- new_r`, works. – David Robinson Aug 02 '12 at 19:51
  • It's a great solution, and probably more elegant than my original code. However, it still doesn't answer the question of where `new_r` was lost. If you get the chance, read my comment on Justin's answer - I'd be interested to get to the bottom of this issue. – MikeZ Aug 03 '12 at 13:25

1 Answers1

6

You're doing a lot of implicit stuff in there... I would suggest writing things more explicitly.

Specifically, mydata.Gest$r <- new_r then replace new_r with r in your plot formula, plot(..., cbind(rs, theo) ~ r, ...). That works for me. Not sure why it works outside of a function and not inside, but relying on plot to look outside the local scope of mydata.Gest for new_r is risky.

Also, use = to assign things to columns in a data frame rather than <-

From a clean session:

data.frame(x<-1:10, y<- 1:10)
ls()

versus

data.frame(x=1:10, y=1:10)
ls()
Justin
  • 42,475
  • 9
  • 93
  • 111
  • Very nice solution, Justin, thank you! I'm still put off by the fact that we haven't pin-pointed where the error is occurring, however. I've used the process I posted before, and I've never had any problems with it. The difference between my previous uses and the one I've posted is that none of the others involve `nls`. The only possible fault I see occurring is that, after using `new_r` to create `GvsR_dataframe`, perhaps the `new_r` object is deleted after being used by `nls`. I am very curious, and would like to get to the bottom of this issue. – MikeZ Aug 03 '12 at 13:18
  • 1
    I think it has more to do with the way `plot.Gest` is written. Its in that call that you run into problems. you can use `new_r` in other portions of the plot call, but the formula portion doesn't seem to look beyond the local scope of the data argument. – Justin Aug 03 '12 at 13:55
  • I've used `plot.Gest` before in a context where the independent variable used was an object outside of the plot function, AND not affiliated with `mydata.Gest`. That's what makes me think it has something to do with `nls`. – MikeZ Aug 03 '12 at 14:14
  • you can add lines to your function like `print(head(new_r))` to help sort this out. you can also use `debugonce(ModelFit)` to step through the function inside its scope. – Justin Aug 03 '12 at 14:21