[R] indexing in a function doesn't work?

[Joshua Wiley]

Hi Benjamin,

You seem to have the right basic ideas, but a lot of your code had
typos and some logic flaws that I guess came from trying to move from
just code to in a function.  I attached the changes I made.  What I
would strongly encourage you to do, is work through each of the little
functions I made and:

1) make sure you understand what it is doing
2) make sure each small function works properly---this means creating
a *variety* of plausible test cases and trying them out
3) once all of the pieces work, then try to wrap them up in your
overall plotter()

The original function you had was large and had many errors, but it is
difficult to debug something when there can be errors coming from
multiple places.  Easier is to break your work into small, tractable
chunks.  Plan in advance what your final goal is, and how each piece
will fit in to that.  Then write each piece and ensure that it works.
>From this point, you will have a much easier time bundling all the
pieces together (even still, there may be additional work, but it will
be more doable because you can be reasonably certain all your code
works, it just does not quite work together.

A few functions I used may be new to you, so I would also suggest
reading the documentation (I know this can be tedious, but it is
valuable)

?match.arg
?switch
?on.exit

note that "..." are arguments passed down from the final function
plotter to internal ones.  They must be named if they are passed to
"...".

You are off to a good start, and I think with some more work, you can
get this going fine.  Long run you may have less headaches and stress
if you take more time at the beginning to write clean code.

I hope this helps,

Josh

On Sun, Apr 1, 2012 at 4:34 PM, Benjamin Caldwell
<btcaldwell at berkeley.edu> wrote:
> Josh,
>
> Many thanks - here's a subset of the data and a couple examples:
>
> plotter(10,3,fram=rwb,framvec=rwb$prcnt.char.depth,obj=prcnt.char.depth,form1=
> post.f.crwn.length~shigo.av,form2=post.f.crwn.length~shigo.av-1,
> form3=leaf.area~(1/exp(shigo.av*x))*n,type=2,xlm=70,ylm=35)
>
> plotter(10,3,fram=rwb, framvec=rwb$prcnt.char.depth, obj=prcnt.char.depth,
> form1= post.f.crwn.length~leaf.area, form2=post.f.crwn.length~leaf.area-1,
> form3=leaf.area~(1/exp(shigo.av*x))*n,type=1, xlm=1500, ylm=35,
> sx=.01,sn=25)
>
>
>
>
> plotter<-function(a,b,fram,framvec,obj,form1,form2,form3, type=1, xlm, ylm,
> sx=.01,sn=25){
> g<-ceiling(a/b)
> par(mfrow=c(b,g))
> num<-rep(0,a)
> sub.plotter<-function(i,fram,framvec,obj,form1,form2,form3,type,
> xlm,ylm,var1,var2){
> temp.i<-fram[framvec <=(i*.10),] #trees in the list that have an attribute
> less than or equal to a progressively larger percentage
> plot(form1, data=temp.i, xlim=c(0,xlm), ylim=c(0,ylm), main=((i-1)*.10))
> if(type==1){
> mod<-lm(form2,data=temp.i)
> r2<-summary(mod)$adj.r.squared
> num[i]<-r2
> legend("bottomright", legend=signif(r2), col="black")
> abline(mod)
> num}
> else{
> if(type==2){
> try(mod<-nls(form3, data=temp.i, start=list(x=sx,n=sn),
> na.action="na.omit"), silent=TRUE)
> try(x1<-summary(mod)$coefficients[1,1], silent=TRUE)
> try(n1<-summary(mod)$coefficients[2,1], silent=TRUE)
> try(lines((1/exp(c(0:70)*x1)*n1)), silent=TRUE)
> try(num[i]<-AIC(mod), silent=TRUE)
> try(legend("bottomright", legend=round(num[i],3) , col="black"),
> silent=TRUE)
> try((num), silent=TRUE)
>   }
> }}
> for(i in 0:a+1){
>  num<-sub.plotter(i,fram,framvec,obj,form1,form2,form3,type,xlm,ylm)
> }
> plot.cor<-function(x){
> temp<-a+1
> lengthx<-c(1:temp)
> plot(x~c(1:temp))
> m2<-lm(x~c(1:temp))
> abline(m2)
> n<-summary(m2)$adj.r.squared
> legend("bottomright", legend=signif(n), col="black")
> slope<-(coef(m2)[2])# slope
> values<-(num)#values for aic or adj r2
> r2ofr2<-(n) #r2 of r2 or AIC
> output<-data.frame(lengthx,slope,values,r2ofr2)
> }
> plot.cor(num)
> write.csv(plot.cor(num)$output,"output.csv") # can't seem to use
> paste(substitute(form3),".csv",sep="") to name it at the moment
> par(mfrow=c(1,1))
> }
>
>
>
>
> On Sun, Apr 1, 2012 at 3:25 PM, Joshua Wiley <jwiley.psych at gmail.com> wrote:
>>
>> Hi,
>>
>> Glancing through your code it was not immediately obvious to me why it
>> does not work, but I can see a lot of things that could be simplified.
>>  It would really help if you could give us a reproducible example.
>> Find/upload/create (in R) some data, and examples of how you would use
>> the function.  Right now, I can only guess what your data etc. are
>> like and based on your description plus what the code you wrote seems
>> to expect to be given.  I could try to give code suggestions, but I
>> have no easy way of testing them so it would be very easy to make
>> typos, etc.  Then you just get back my edits to your code that still
>> do not work and maybe it is because of something fundamentally wrong
>> with what I have done, a simple typo, or something else still wrong in
>> your code that I did not fix.
>>
>> Anyway, if you send some data and an example using your function
>> (i.e., using the data you send, write our form1, form2, type, etc.), I
>> will take a look at your function and see if I can make it run.
>>
>> Cheers,
>>
>> Josh
>>
>> On Sun, Apr 1, 2012 at 3:13 PM, Benjamin Caldwell
>> <btcaldwell at berkeley.edu> wrote:
>> > Hello,
>> >
>> > I've written a small function that's supposed to save me some time, and
>> > it's ending up killing it- the intention is to iteratively subset a
>> > dataset
>> > fram on framevec, fit a model (either lm or nls depending on type) and
>> > return the r2 or AIC from the model, respectively. Although as far as I
>> > can
>> > tell in my code the plots are dependent on the fit of the model to the
>> > data
>> > and the r2 and AIC reported are also dependent on the re-fitted model,
>> > the
>> > plots show the same linear or non-linear model used for every subset of
>> > the
>> > data.
>> >
>> > However, the r2 and AIC values come back different for each subset.
>> >
>> > When I do the subsetting and model fitting outside the function, the
>> > model
>> > fit is different, with different slopes, for each subset of the data.
>> >
>> > I'm going to end up doing this without the function if I don't solve
>> > this
>> > soon. Any help much appreciated.
>> >
>> >
>> > #a is the number of times to loop the tenth of percent, b is first
>> > dimension of mfrow, frame is the dataframe, framevec is the vector
>> > within
>> > the dataframe you're using to subset (should be a percentage), form 1 is
>> > the simple y~x for the plot, form 2 is y~x for regression, type is lm 1
>> > or
>> > 2 nls ,form 3 is the formula for the nls, sx and sn are the start values
>> > for nls
>> >
>> > plotter<-function(a,b,fram,framvec,form1,form2,form3, type=1, xlm, ylm,
>> > sx=.01,sn=25){
>> > g<-ceiling(a/b)
>> > par(mfrow=c(b,g))
>> >  num<-rep(0,a)
>> > sub.plotter<-function(i,fram,framvec,form1,form2,form3,type,
>> > xlm,ylm,var1,var2){
>> > temp.i<-fram[framvec <=(i*.10),]
>> >  plot(form1, data=temp.i, xlim=c(0,xlm), ylim=c(0,ylm),
>> > main=((i-1)*.10))
>> > if(type==1){
>> >  mod<-lm(form2,data=temp.i)
>> > r2<-summary(mod)$adj.r.squared
>> > num<-r2
>> >  legend("bottomright", legend=signif(r2), col="black")
>> > abline(mod)
>> >  num}
>> > else{
>> > if(type==2){
>> >  try(mod<-nls(form3, data=temp.i, start=list(x=sx,n=sn),
>> > na.action="na.omit"), silent=TRUE)
>> > try(x1<-summary(mod)$coefficients[1,1], silent=TRUE)
>> >  try(n1<-summary(mod)$coefficients[2,1], silent=TRUE)
>> > try(lines((1/exp(c(0:70)*x1)*n1)), silent=TRUE)
>> >  try(num[i]<-AIC(mod), silent=TRUE)
>> > try(legend("bottomright", legend=round(num[i],3) , col="black"),
>> > silent=TRUE)
>> >  try((num), silent=TRUE)
>> >  }
>> > }}
>> > for(i in 0:a+1){
>> >  num<-sub.plotter(i,fram,framvec,form1,form2,form3,type,xlm,ylm)
>> > }
>> > plot.cor<-function(x){
>> > temp<-a+1
>> > lengthx<-c(1:temp)
>> >  plot(x~c(1:temp))
>> > m2<-lm(x~c(1:temp))
>> > abline(m2)
>> >  n<-summary(m2)$adj.r.squared
>> > legend("bottomright", legend=signif(n), col="black")
>> >  slope<-(coef(m2)[2])# slope
>> > values<-(num)#values for aic or adj r2
>> > r2ofr2<-(n) #r2 of r2 or AIC
>> >  output<-data.frame(lengthx,slope,values,r2ofr2)
>> > }
>> > plot.cor(num)
>> > write.csv(plot.cor(num)$output,"output.csv") # can't seem to use
>> > paste(substitute(form3),".csv",sep="") to name it at the moment
>> > par(mfrow=c(1,1))
>> > }
>> >
>> > Ben
>> >
>> >        [[alternative HTML version deleted]]
>> >
>> > ______________________________________________
>> > R-help at r-project.org mailing list
>> > https://stat.ethz.ch/mailman/listinfo/r-help
>> > PLEASE do read the posting guide
>> > http://www.R-project.org/posting-guide.html
>> > and provide commented, minimal, self-contained, reproducible code.
>>
>>
>>
>> --
>> Joshua Wiley
>> Ph.D. Student, Health Psychology
>> Programmer Analyst II, Statistical Consulting Group
>> University of California, Los Angeles
>> https://joshuawiley.com/
>
>



-- 
Joshua Wiley
Ph.D. Student, Health Psychology
Programmer Analyst II, Statistical Consulting Group
University of California, Los Angeles
https://joshuawiley.com/
-------------- next part --------------
## your example csv
rwb <- read.csv(file.choose())

## linear model function
bclm <- function(f, data) {
  mod <- lm(f, data=data)
  r2 <- summary(mod)$adj.r.squared
  legend("bottomright", legend=signif(r2), col="black")
  abline(mod)
  return(r2)
}

## nonlinear model function
bcnls <- function(f, data, start) {
  mod <- nls(f, data=data, start=start, na.action="na.omit")
  x1 <- summary(mod)$coefficients[1,1]
  n1 <- summary(mod)$coefficients[2,1]
  lines((1/exp((0:70)*x1)*n1))
  aic <- AIC(mod)
  legend("bottomright", legend=round(aic,3) , col="black")
  return(aic)
}

## subplotter creates a plot and depending on whether type = "lm" or type = "nls"
## uses the linear model function (bclm) or nonlinear (bcnls)
sub.plotter <- function(i, fram, framvec, obj, form1, form2, form3, type = c("lm", "nls"), xlm, ylm, s = list(x=.01, n=25)) {
  # trees in the list that have an attribute less than or equal to a progressively larger percentage
  temp.i <- fram[framvec <= (i*.10), ]
  plot(form1, data=temp.i, xlim = c(0, xlm), ylim = c(0, ylm), main = ((i-1)*.10))

  type <- match.arg(type)

  out <- switch(type,
    lm = bclm(form1, temp.i),
    nls = try(bcnls(form3, temp.i, start = s), silent = TRUE))
  if (!is.numeric(out)) {out <- NA}
  return(out)
}

plot.cor <- function(x, a) {
  lengthx <- seq.int(a + 1)
  plot(x ~ lengthx)
  m2 <- lm(x ~ lengthx)
  abline(m2)
  n <- summary(m2)$adj.r.squared
  legend("bottomright", legend = signif(n), col = "black")
  slope <- coef(m2)[2]# slope
  values <- x #values for aic or adj r2
  r2ofr2 <- n #r2 of r2 or AIC
  output <- data.frame(lengthx, slope, values, r2ofr2)
  return(output)
}

plotter <- function(a, b, ...) {
  oldpar <- par(no.readonly = TRUE)
  on.exit(par(oldpar))

  g <- ceiling(a/b)
  par(mfrow = c(b, g))
  num <- vector("numeric", a)
  for(i in (0:a)+1) {
    num[i] <- sub.plotter(i, ...)
  }
  output <- plot.cor(num, a)
  write.csv(output, "output.csv")
}

## final function
plotter(a = 10, b = 3,
  fram=rwb, framvec=rwb$prcnt.char.depth, obj=rwb$prcnt.char.depth,
  form1 = post.f.crwn.length ~ shigo.av,
  form2 = post.f.crwn.length ~ shigo.av-1,
  form3 = leaf.area ~ (1/exp(shigo.av*x))*n,
  type="nls", xlm=70, ylm=35)