#a modification to the WIC, BIC and TNW calculations of Roughgarden (1972, Am Nat) and Bolnick et al (2002, Ecology), and in the IndSpec program. These calculations work for continuous diet data such as prey sizes when the number of diet items varies among individuals. 
#Reference: Ingram, T., W.E. Stutz & D.I. Bolnick. 2011. Does intraspecific size variation in a predator affects its diet diversity and top-down control of prey? PLoS ONE 6: e20782.

#file format: input a matrix where each row is an individual with one prey size per column, and NA's filling up the rest of the matrix. The number of columns should be the maximum number of items consumed by any individual.
#example:
#1.5,0.5, NA, NA, NA
#0.5,1.0,2.5,3.0,1.5
#2.0,1.5,0.5, NA, NA

#simulating data
Nind<-10 #number of individuals
mean_Nitems<-5 #average # of prey items/individual
approxIS<-0.5 #tune the degree of diet variation among individuals - may not correspond to the true expected simulated values, but illustrates that the calculations work in populations with high or low individual specialization. 
Ns<-rpois(Nind,lambda=mean_Nitems-1)+1 # +1 to ensure no zeros 
d<-matrix(NA,nrow=Nind,ncol=max(Ns))
for(i in 1:Nind)d[i,1:Ns[i]]<-rnorm(1,0,sqrt(1-approxIS))+rnorm(Ns[i],0,sqrt(approxIS))

#first read in the functions 'IndSpecContinuous' and 'Var' below

IS_eq<-IndSpecContinuous(diet=d,weights="equal");IS_eq
IS_Ni<-IndSpecContinuous(diet=d,weights="N_items");IS_Ni

#confirm that BIC+WIC=TNW (might not be precisely equal depending on machine precision)
sum(IS_eq[1:2]);IS_eq[3]
sum(IS_Ni[1:2]);IS_Ni[3]

#main function
IndSpecContinuous<-function(diet,weights="equal"){

if(weights%in%c("equal","N_items")==FALSE)stop("`weights` must be either `equal` or `N_items`")

#number of items consumed by each individual
Ns<-apply(diet,1,function(x)sum(!is.na(x)))
indvars<-apply(diet,1,Var)
indmeans<-apply(diet,1,mean,na.rm=T)
meansize<-mean(indmeans)

#weighting by the number of items in each individual's diet, so those with more data contribute more to estimating parameters
if(weights=="N_items"){
	BIC<-Var(indmeans,w=Ns)
	WIC<-weighted.mean(indvars,w=Ns,na.rm=T)
	TNW<-Var(as.numeric(diet))
}else{
#weighting each fish equally regardless of diet item number; requires weighting each diet item by the inverse of the number of items in the individual's diet for calculating TNW
	wts<-diet
	for(i in 1:dim(diet)[1])wts[i,1:Ns[i]]<-1/Ns[i]
	BIC<-Var(indmeans)
	WIC<-mean(indvars)
	TNW<-Var(as.numeric(diet),w=as.vector(wts),w.mean=TRUE)
	}
c(BIC=BIC,WIC=WIC,TNW=TNW,WIC_TNW=WIC/TNW)
}

#helper function - calculates (optionally weighted) maximum likelihood variances (not the unbiased estimator)
Var<-function(x,w=rep(1/length(x),length(x)),w.mean=TRUE){
	w<-w[!is.na(x)];x<-x[!is.na(x)]
	if(length(x)==0){
		result<-NA
		}else{
	if(sum(w)!=1)w<-w/sum(w)
	meanx<-ifelse(w.mean==TRUE,weighted.mean(x,w=w),mean(x))
	result<-sum(((x-meanx)^2)*w)
	}	}