Plotting problem in R

#1
Hi,

I am relatively new to R and have encountered a plotting problem last night that I am unable to solve. Before last night, this same exact code worked, but since then I've tried it on multiple computers with no success. I am using EPA-published code to build a predictive model for my thesis. The error produced as well as the code are provided below.

Error in plot.window(...) : need finite 'xlim' values
In addition: Warning messages:
1: In min(x) : no non-missing arguments to min; returning Inf
2: In max(x) : no non-missing arguments to max; returning -Inf


predall<-read.table("WSU_predictors_working_HUC12_Arc_Sine_transformation_10.24.txt",row.names="Sample",header=T,sep="\t");
predall<-read.table("WSU_predictors_working_HUC12_Arc_Sine_transformation_10.24.txt",row.names="Sample_ID",header=T,sep="\t");
head(predall);
dim(predall);
candvar<-c("HUC12_AREA_HA_log10","ELEV_m","M_Slp_sqt","Precip_mm","Temp_CX10","Lith_R.1._E.0.","Baseflow","Strm_Power_log10","Road_Density","Forest_Fragmentation","LC_Diversity_Site","Linear_Dist_Dam_km","Linear_Dist_NPDES_km","Linear_Dist_Fish_Barrier_km","ASIN_Percent_Urban_HUC12","ASIN_Percent_Urban_Per_50m_Buffer_HUC12","ASIN_Percent_Forest_HUC12","ASIN_Percent_Forest_Per_50m_Buffer_HUC12","ASIN_Percent_Shrub_Scrub_HUC12","ASIN_Percent_Shrub_Scrub_Per_50m_Buffer_HUC12","ASIN_Percent_Agric_HUC12","ASIN_Percent_Agric_Per_50m_Buffer_HUC12","ASIN_Percent_Wetl_HUC12","ASIN_Percent_Wetl_Per_50m_Buffer_HUC12");
windows();
sapply(candvar,function(nam)hist(predall[,nam],main=nam));
bugall<-read.table("WSU_sitspp-tab.txt",row.names="Sample_ID",header=T,sep="\t");
bugall<-bugall[row.names(predall),];
row.names(bugall)==row.names(predall);
bugall.pa<-bugall;
bugall.pa[bugall.pa>0]<-1;
predcal<-predall[predall[,'Reference_Test']=='R',];
pred.vld<-predall[substr(as.character(predall[,'Reference_Test']),1,1)=='R',];
bugcal<-bugall[predall[,'Reference_Test']=='R',];
bugcal.pa<-bugall.pa[predall[,'Reference_Test']=='R',];
bug.vld.pa<-bugall.pa[substr(as.character(predall[,'Reference_Test']),1,1)=='R',];
psite.occ<-colSums(bugcal.pa)/dim(bugcal.pa)[[1]];
nonrare.taxa<-names(psite.occ)[psite.occ>0.05];
bugcal.nonrare<- bugcal[,nonrare.taxa];
bugcal.pa.nonrare<- bugcal.pa[,nonrare.taxa];
source("dapply.r");
sornfun<-function(siti,sitj) {
shared<-sum((siti>0)&(sitj>0));
uniquei<-sum((siti>0)&(sitj==0));
uniquej<-sum((siti==0)&(sitj>0));
1-(2*shared/(2*shared+uniquei+uniquej)); #return Sorenson dissimilarity;
} #end of function;
dissim<-dapply(bugcal.pa.nonrare,1,bugcal.pa.nonrare,1,sornfun);
clus1<-agnes(x=dissim,diss=T,method="flexible", par.method=0.8,keep.diss=F,keep.data=F);
plot(clus1,which.plots=2,labels=row.names(bugcal.pa),cex=.4);
ccc<-identify(as.hclust(clus1)); #interactive ID of clusters on dendrogram;
grp.5<-rep(NA,dim(bugcal)[[1]]);
for(k in 1:length(ccc)) {grp.5[ccc[[k]]]<-k};
table(grp.5); #count number of sites in each group;
cbind(row.names(predcal),grp.5); #list calibration sites and their group assignments;
nkeep<-c(rep(5,6),1);
source("dfa.allsub.v4.r");
start.time=proc.time();
dfm.best<-dfa.allsub.v3(bug.cal=bugcal.pa,bug.vld=bug.vld.pa,pred.cal=predcal,pred.vld=pred.vld,grps=grp.5,candvar=candvar,numkeep=nkeep,Pc=0.5);
elaps<-proc.time()-start.time;
print(c("elapsed time = ",elaps));
dfm.best$null.stats;
bestmods<-dfm.best
bestmods
par(mfrow=c(1,1));
plot(bestmods$order,bestmods$RMSE.cal,ylim=c(0.15,0.20),type='p',pch='C', col='blue',
cex=.7,xlab='Model order',ylab='RMSE(O/E)');
par(mfrow=c(1,1));
plot(bestmods$order,bestmods$RMSE.cal,ylim=c(1,7),ylim=c(0.15,0.20),type='p',pch='C', col='blue',
cex=.7,xlab='Model order',ylab='RMSE(O/E)');
plot(bestmods$order,bestmods$cls.crct.resub,ylim=c(20,60),type='p',pch='R',
cex=.5,xlab='Model order',ylab='Percent correct');
plot(bestmods$order,bestmods$RMSE.cal,xlim=c(1,7),ylim=c(0.15,0.20),type='p',pch='C', col='blue',
cex=.5,xlab='Model order',ylab='Percent correct');
par(mfrow=c(1,1));
plot(bestmods$order,bestmods$RMSE.cal,xlim=c(1,7),ylim=c(0.15,0.20),type='p',pch='C', col='blue',
cex=.5,xlab='Model order',ylab='Percent correct');
points(bestmods$order,bestmods$cls.crct.cv,pch='C',cex=.5);
legend(locator(1),legend=c('Resubstitution','Crossvalidation'),pch=c('R','C'));
plot(bestmods$order,bestmods$cls.crct.resub,xlim=c(1,7),ylim=c(20,60),type='p',pch='R',
cex=.5,xlab='Model order',ylab='Percent correct');
points(bestmods$order,bestmods$cls.crct.cv,pch='C',cex=.5);
legend(locator(1),legend=c('Resubstitution','Crossvalidation'),pch=c('R','C'));
par(mfrow=c(1,1));
plot(bestmods$order,bestmods$RMSE.cal,xlim=c(1,7),ylim=c(0.15,0.20),type='p',pch='C', col='blue',
cex=.7,xlab='Model order',ylab='RMSE(O/E)');
points(bestmods$order,bestmods$RMSE.vld,pch='V',cex=.7,col='red');
abline(dfm.best$null.stats["RMSE.cal"],0,lty=1);
abline(dfm.best$null.stats["RMSE.vld"],0,lty=2);


As you can see, I tried adding the xlim value (which was not previously used and the plot worked just fine). but all this produced is a blank graph with labels. I appreciate any help I can get.

Thanks in advance!
 

bugman

Super Moderator
#2
I deleted your other post - please just post in one forum, it makes it easier to keep track on replies and ensures ansewers are not repeated by different members.