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/
>
______________________________________________
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.
