Josh
Apologies I haven't responded earlier. This looks great - I ended up doing
what I needed done piece-by-piece because of a looming deadline, but
understanding this code and your suggestions below will be a weekend
project.
Many thanks for all your help.
*
*
Best
Ben
On Tue, Apr 3, 2012 at 8:01 AM, Joshua Wiley <jwiley.ps...@gmail.com> wrote:
> 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
> <btcaldw...@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.ps...@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
> >> <btcaldw...@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@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/
>
[[alternative HTML version deleted]]
______________________________________________
R-help@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.
