EBSeq/R/PolyFitPlot.R

   1 PolyFitPlot <-
   2 function(X , Y , nterms , xname="Estimated Mean", yname="Estimated Var", pdfname="", xlim=c(-1,5), ylim=c(-1,7), ChangeXY=F,col="red"){
   3
   4         b=rep(NA,nterms)
   5         logX=matrix(rep(X, nterms),ncol=nterms, byrow=T)
   6         for (i in 1:nterms)
   7                 logX[,i]=(log10(X))^i
   8         colnames(logX)=paste("logmu^",c(1:nterms))
   9         rownames(logX)=names(X)
  10         NotUse=c(names(X)[X==0],names(Y)[Y==0],names(X)[rowMeans(logX)==-Inf],names(X)[rowMeans(logX)==Inf])
  11         Use=names(X[!names(X)%in%NotUse])
  12         Lm=lm(log10(Y[Use])~logX[Use,1:nterms])
  13         b=summary(Lm)$coefficients[2:(nterms+1),1]
  14         d=summary(Lm)$coefficients[1,1]
  15         bvec=matrix(rep(b,length(X)),ncol=nterms,byrow=T)
  16         fit=rowSums(logX*bvec)+d
  17         main2=NULL
  18         if (ChangeXY==T){
  19                 X.plot=log10(Y)
  20                 Y.plot=log10(X)
  21                 fit.X.plot=fit
  22                 fit.Y.plot=log10(X)
  23         }
  24         else{
  25         X.plot=log10(X)
  26         Y.plot=log10(Y)
  27             fit.X.plot=log10(X)
  28                 fit.Y.plot=fit
  29                                    }
  30
  31         for (i in 1:nterms)
  32                 main2=paste(main2,round(b[i],2),"*log(",xname,")^",i,"+")
  33         main=pdfname
  34
  35         smoothScatter(X.plot, Y.plot ,main=main,xlim=xlim,ylim=ylim,xlab=xname,ylab=yname,axes=F)
  36         axis(1,at=seq(xlim[1],xlim[2],by=1), 10^seq(xlim[1],xlim[2],by=1))
  37         axis(2,at=seq(ylim[1],ylim[2],by=2), 10^seq(ylim[1],ylim[2],by=2))
  38         Sortit=order(fit.X.plot)
  39         lines(fit.X.plot[Sortit],fit.Y.plot[Sortit],col=col,lwd=3)
  40         output=list(b=b,d=d,lm=Lm,fit=fit,sort=Sortit)
  41         names(output$b)=paste(xname,"^",c(1:length(output$b)))
  42         output
  43 }
  44