mc_int

mc_int.r

mc.int<-function(mod.in,int=NULL,n=10000,int.base=0,plot=T,plot.pts=T,
  plot.poly=T,cols=list('black','green','blue'),round.val=2,val=F){
	
	require(scales)
	
	#variable names from model
	if(is.null(mod.in$call$formula)){
		x.names<-colnames(eval(mod.in$call$x))
		y.names<-colnames(eval(mod.in$call$y))
		x.range<-range(eval(mod.in$call$x))
		}
	else{
		forms<-eval(mod.in$call$formula)
		dat.names<-model.frame(forms,data=eval(mod.in$call$data))
		y.names<-as.character(forms)[2]
		x.names<-names(dat.names)[!names(dat.names) %in% y.names]
		x.range<-range(eval(mod.in$call$data)[,x.names])
		}
	
	#integration range if null
	if(is.null(int)) int<-x.range
	
	#create data for model fit
	x.fit<-seq(x.range[1],x.range[2],length=n)
	pred.dat<-data.frame(x.fit)
	names(pred.dat)<-x.names
	y.fit<-predict(mod.in, newdata=pred.dat)
	dat.fit<-data.frame(x.fit,y.fit)
	
	#create random integration data
	x.rand<-runif(n,int[1],int[2])
	y.rand<-runif(n,int.base,max(y.fit))
	dat.rand<-data.frame(x.rand,y.rand,y.under=rep(cols[[1]],n),stringsAsFactors=F)
	names(dat.rand)[1]<-x.names
	dat.rand$y.under[y.rand<predict(mod.in,newdata=dat.rand[,1,drop=F])]<-cols[[2]]
	y.under<-dat.rand$y.under
	
	#get integration
	area.all<-diff(range(x.rand))*diff(range(y.rand))
	area.und<-area.all*sum(y.under==cols[[2]])/length(y.under)
	out.txt<-paste('Integration from',int[1],'to',int[2],'=',round(area.und,round.val),
		sep=' ')
	
	if(plot){
		
		plot(x.rand,y.rand,xlim=x.range,ylim=c(int.base,max(y.fit)),type='n',main=out.txt)
		points(dat.fit$x.fit,dat.fit$y.fit,type='l')

		if(plot.poly){
			
			#create data for integration polygon
			cord.x<-seq(int[1],int[2],length=(n*diff(int)/diff(x.range)))
			cord.x<-data.frame(cord.x)
			names(cord.x)<-x.names
			cord.y<-c(int.base,predict(model,newdata=cord.x),int.base)
			cord.x<-rbind(int[1],cord.x,int[2])
			dat.poly<-data.frame(cord.x,cord.y)
			names(dat.poly)[1]<-x.names
			
			polygon(dat.poly[,1],dat.poly[,2],col=alpha(cols[[3]],0.5),border=NA)
			
			}
		
		if(plot.pts){
			
			points(dat.rand[,1],dat.rand[,2],col=alpha(dat.rand$y.under,0.5),pch=19)
			
			}

		}

	if(!val) return(cat(out.txt))
	area.und

	}