R/alkD.R

In fishmethods: Fishery Science Methods and Models

alkD<-function(x,lss=NULL, minss=NULL,maxss=NULL,sampint=NULL,allocate=1){
      if(is.null(x)) 
         stop ("No data")
      if(is.null(lss)) 
         stop ("length sample size is required") 
      if(is.null(minss)) 
         stop ("minimum age is required") 
      if(is.null(maxss)) 
         stop ("maximum age is required") 
      if(is.null(sampint)) 
         stop ("sample size interval is required") 
      if(maxss>lss) 
         stop ("maxss can not be larger than lss") 

      nages<-ncol(x)-2
	alpha<-x[,2]/sum(x[,2])				
	x$sumA<-rowSums(x[,3:as.numeric(nages+2)])
	thetala<-x                                     
	for(i in 1:length(thetala[,1])){                                            
        for(j in 3:as.numeric(nages+2)){
           if(thetala$sumA[i]>0) thetala[i,j]<-thetala[i,j]/thetala$sumA[i]
           if(thetala$sumA[i]==0) thetala[i,j]<-0

        }
      }

	theta<-cbind(thetala,alpha)
	ta<-theta
	for(i in 1:length(ta[,1])){                                             
  	   for(j in 3:as.numeric(nages+2)){
            ta[i,j]<-ta[i,j]*ta$alpha[i]
         }
      }
	est_theta_a<-colSums(ta[,3:as.numeric(nages+2)])

	if(allocate==1){
 	  var1<-theta
        for(i in 1:length(var1[,1])){                                            
    		for(j in 3:as.numeric(nages+2)){
                var1[i,j]<-var1[i,j]*(1-var1[i,j])*var1$alpha[i]
            }
        }
       Va<-colSums(var1[,3:as.numeric(nages+2)])
     }
	if(allocate==2){
 	   var1<-theta
  	   for(i in 1:length(var1[,1])){                                            
            for(j in 3:as.numeric(nages+2)){
                var1[i,j]<-var1[i,j]*(1-var1[i,j])*var1$alpha[i]^2
            }
         }
        Va<-colSums(var1[,3:as.numeric(nages+2)])*length(var1[,2])
      }

	var2<-theta
	for(j in 3:as.numeric(nages+2)){                                        
         for(i in 1:length(var2[,1])){
            var2[i,j]<-var2$alpha[i]*(var2[i,j]-est_theta_a[j-2])^2
         }
         Ba<-colSums(var2[,3:as.numeric(nages+2)])

      }
	ageSS<-seq(minss,maxss,sampint)               
	D<-round(sqrt(sum(Va)/ageSS+sum(Ba)/lss),4)	
      dd<-as.data.frame(cbind(ageSS,D))
      names(dd)<-c("age sample size","D")
      SS<-c(paste("D statistic for L=",lss,ifelse(allocate==1,
           ", Allocation: Proportional",", Allocation: Fixed")))
      output<-list(SS,dd);names(output)<-c("label","results")
      return(output)	
 }