Condition_indicator.Rmd

# Fish Condition Indicator

## Contributor name(s)
Laurel Smith

## Data name
Relative condition

##Methods
Relative condition was introduced by @Cren1951a as a way to remove the influence of length on condition, and @Blackwell2000 noted that Kn may be useful in detecting prolonged physical stress on a fish populations. Relative condition is calculated as:
Kn = W/W'
Where W is the weight of an individual fish and W' is the predicted length-specific mean weight for the fish population in a given region. 
Here, relative condition was calculated for finfish stocks commonly caught on the Northeast Fisheries Science Center's (NEFSC) autumn and spring bottom trawl surveys. Where data allowed, length-weight parameters were calculated for W' by species, sex and season over the time period 1992-20XX. When sample sizes of individual fish weights and lengths were too low, parameters were calculated for aggregated spring and fall survey data over the same time period. Fall survey relative condition was calculated by sex for those species that exhibited differences in growth between sexes and aggregated across sex for those that did not.

### Data source(s)
Individual fish lengths (nearest 0.5 cm) and weights (grams) were collected on the NEFSC bottom trawl surveys from 1992-20XX aboard RVs Albatross IV, Deleware II and the Henry B. Bigelow. A small number of outlier values were removed when calculating the length-weight parameters, with nominal effect on the results.

### Data extraction
Data were extracted from NEFC's survey database (SVDBS) using SQL. 

Write SQL query here
```{sql, eval = F, echo = T}
SELECT cruise6,stratum,tow,station,
  year,month,day,time,beglat,beglon,setdepth,
  surftemp,bottemp,
  svspp,sex,length,age,maturity,indid,indwt,stom_volume,stom_wgt, expcatchwt, expcatchnum
from connection to oracle
(select b.cruise6,b.stratum,b.tow,b.station,
  s.est_year year,est_month month,est_day day,
  substr(est_time,1,2)||substr(est_time,4,2) time,
  round(substr(beglat,1,2) + (substr(beglat,3,7)/60),6) beglat,
  round(((substr(beglon,1,2) + (substr(beglon,3,7)/60)) * -1), 6) beglon,
  setdepth,surftemp, bottemp,
  b.svspp,sex,length,age,maturity,indid,indwt,stom_volume,stom_wgt, expcatchwt, expcatchnum
from union_fscs_svbio b, union_fscs_svcat p, union_fscs_svsta s, svdbs_cruises c
where 
  season = &sson and
    b.svspp in ('013','015','023','026','028','032','072','073','074','075','076','077','078','102','103','104','105','106','107','108','121','131','135','141','143','145','155','164','193','197') and
  (b.cruise6=s.cruise6) and
  (c.cruise6=b.cruise6) and
  (p.cruise6=c.cruise6) and
  (p.stratum=b.stratum) and
  (b.stratum=s.stratum) and
  (p.station=b.station) and
  (b.station=s.station) and
  (p.svspp=b.svspp) and
  (p.tow=b.tow) and
  (b.tow=s.tow) );

  %put &sqlxmsg;
  %put &sqlxrc;

create view spp as

select comname, svspp
from connection to oracle
(select comname, svspp
from svspecies_list);

  %put &sqlxmsg;
  %put &sqlxrc;

execute (commit) by oracle;

```

### Data analysis
The following growth curve was fit through individual fish lengths and weights from the NEFSC bottom trawl survey data from 1992-20XX to produce the reference length-weight parameters in R:


$$\textrm{Weight} = e^{Fall_{coef}} * \textrm{Length}^{Fall_{exp}}$$

Where length is in cm and weight is in kg. Fall survey data were used where sample sizes allowed for growth curve estimation, otherwise data from spring and fall seasons were combined. 
Individual fish lengths from NEFSC fall bottom trawl survey from 1992-2017 were then used to calculate predicted weights using the reference length-weight parameters. Relative condition (Kn) was calculated annually by species and sex (for sexually dimorphic species) by dividing individual fish weights by the predicted weight. 


Include accompanying R code here:
```{r analysis, echo = T, eval = F}
# Length-weight parameter calculation:
function (data, min.n = 25, min.range = 5, data.avail = NA, data.avail.bigelow = NA) 
{
	if(is.null(dim(data.avail))) data.avail <- lw.data.availability(data, min.n, min.range)
	data.avail <- data.avail[apply(data.avail[,2:5], 1, any),]
	if(is.null(dim(data.avail.bigelow)))data.avail.bigelow <- lw.data.availability(data[data$data.source == "Bigelow",], min.n, min.range)
	data.avail.bigelow <- data.avail.bigelow[apply(data.avail.bigelow[,2:5], 1, any),]
	data.spp <- as.numeric(rownames(data.avail[data.avail$sex.season == TRUE,]))
	lw.output <- data.frame(matrix(ncol = 12))
	names(lw.output) <- c("species.name", "species.code", "source", "sex", "season", "slope", "slope.p", "intercept", "intercept.p", "min.length", "max.length", "check.diff")
	for(sp in data.spp){
		sp.data <- lw.data[lw.data$species == sp,]
		sp.name <- unique(as.character(species.names$scientific_name[species.names$svspp == sp]))
		print(sp.name)
#All model
		print("Species Level")
		this.data <- bigelow.test(sp.data, data.avail[as.character(sp),], data.avail.bigelow[as.character(sp),], "weight.log~length.log", "species")
		ds <- this.data[[2]]
		this.data <- this.data[[1]]
		if(!is.null(dim(this.data))){
			this.model <- lm(weight.log~length.log, data = this.data)
			model.coefs <- coef(this.model)
			model.summary <- coef(summary(this.model))
			length.range <- range(sp.data$length)
			length.log <- log(seq(length.range[1], length.range[2], by=.5))
			species <- rep(sp, length(length.log))
			predict.length <- data.frame(species, length.log)
			check.diffs <- plot.lw(this.data, this.model, "species", sp.name, predict.length)
			lw.output <- rbind(lw.output, c(sp.name, sp, ds, NA, NA, model.coefs["length.log"], model.summary["length.log", "Pr(>|t|)"], model.coefs["(Intercept)"], model.summary["(Intercept)", "Pr(>|t|)"], length.range[1], length.range[2], check.diffs))
		}
#Sex model
		print("Sex Level")
		this.data <- sp.data[sp.data$sex > 0,]
		model.definition <- "weight.log~length.log * factor(sex)"
		this.data <- bigelow.test(this.data, data.avail[as.character(sp),], data.avail.bigelow[as.character(sp),], model.definition, "sex")
		ds <- this.data[[2]]
		this.data <- this.data[[1]]
		if(!is.null(dim(this.data))){
			this.model <- lm(weight.log~length.log * factor(sex), data = this.data)
			male.range <- range(this.data$length[this.data$sex == 1])
			female.range <- range(this.data$length[this.data$sex == 2])			
			model.summary <- coef(summary(this.model))
			if(any(model.summary[grep("sex", rownames(model.summary)), "Pr(>|t|)"] <= .05)){
				model.coefs <- coef(this.model)
				length.log <- rep(log(seq(length.range[1], length.range[2], by=.5)),2)
				sex <- c(rep(1, length(length.log)/2), rep(2, length(length.log)/2))
				predict.length <- data.frame(length.log, sex)
				check.diffs <- plot.lw(this.data, this.model, "sex", sp.name, predict.length)
				lw.output <- rbind(lw.output, c(sp.name, sp, ds, 1, NA, model.coefs["length.log"], model.summary["length.log", "Pr(>|t|)"], model.coefs["(Intercept)"], model.summary["(Intercept)", "Pr(>|t|)"], male.range[1], male.range[2], check.diffs[1]))
				lw.output <- rbind(lw.output, c(sp.name, sp, ds, 2, NA, model.coefs["length.log"] + model.coefs["length.log:factor(sex)2"], model.summary["length.log:factor(sex)2", "Pr(>|t|)"], model.coefs["(Intercept)"] + model.coefs["factor(sex)2"], model.summary["factor(sex)2", "Pr(>|t|)"], female.range[1], female.range[2], check.diffs[2]))
			}
			else{
				print(paste("Model parameters not significantly different for", model.definition))
			}
		}
#Season model
		model.definition <- "weight.log~length.log * factor(season)"
		this.data <- bigelow.test(sp.data, data.avail[as.character(sp),], data.avail.bigelow[as.character(sp),], model.definition, "season")
		ds <- this.data[[2]]
		this.data <- this.data[[1]]
		if(!is.null(dim(this.data))){
			this.model <- lm(weight.log~length.log * factor(season), data = this.data)
			fall.range <- range(this.data$length[this.data$season == "FALL"])
			spring.range <- range(this.data$length[this.data$season == "SPRING"])
			model.summary <- coef(summary(this.model))
			if(any(model.summary[grep("season", rownames(model.summary)), "Pr(>|t|)"] <= .05)){
				model.coefs <- coef(this.model)
				length.log <- rep(log(seq(length.range[1], length.range[2], by=.5)),2)
				season <- c(rep("FALL", length(length.log)/2), rep("SPRING", length(length.log)/2))
				predict.length <- data.frame(length.log, season)
				check.diffs <- plot.lw(this.data, this.model, "season", sp.name, predict.length)
				lw.output <- rbind(lw.output, c(sp.name, sp, ds, NA, "FALL", model.coefs["length.log"], model.summary["length.log", "Pr(>|t|)"], model.coefs["(Intercept)"], model.summary["(Intercept)", "Pr(>|t|)"], fall.range[1], fall.range[2], check.diffs[1]))
				lw.output <- rbind(lw.output, c(sp.name, sp, ds, NA, "SPRING", model.coefs["length.log"] + model.coefs["length.log:factor(season)SPRING"], model.summary["length.log:factor(season)SPRING", "Pr(>|t|)"], model.coefs["(Intercept)"] + model.coefs["factor(season)SPRING"], model.summary["factor(season)SPRING", "Pr(>|t|)"], spring.range[1], spring.range[2], check.diffs[1]))
			}
			else{
				print(paste("Model parameters not significantly different for", model.definition))
			}
		}
#Sex-Season model
		this.data <- sp.data[sp.data$sex > 0,]
		model.definition <- "weight.log~length.log * factor(sex) * factor(season)"
		this.data <- bigelow.test(this.data, data.avail[as.character(sp),], data.avail.bigelow[as.character(sp),], model.definition, c("sex","season"))
		ds <- this.data[[2]]
		this.data <- this.data[[1]]
		if(!is.null(dim(this.data))){
			this.model <- lm(weight.log~length.log * factor(sex) * factor(season), data = this.data)
			model.summary <- coef(summary(this.model))
			male.fall.range <- range(this.data$length[this.data$season == "FALL" & this.data$sex == 1])
			male.spring.range <- range(this.data$length[this.data$season == "SPRING" & this.data$sex == 1])
			female.fall.range <- range(this.data$length[this.data$season == "FALL" & this.data$sex == 2])
			female.spring.range <- range(this.data$length[this.data$season == "SPRING" & this.data$sex == 2])
			if(any(model.summary[grep("season", rownames(model.summary)), "Pr(>|t|)"] <= .05 | any(model.summary[grep("season", rownames(model.summary)), "Pr(>|t|)"] <= .05))){
				model.coefs <- coef(this.model)
				male.fall.int <- model.coefs["(Intercept)"]
				male.fall.int.p <- model.summary["(Intercept)", "Pr(>|t|)"]
				male.fall.slope <- model.coefs["length.log"]
				male.fall.slope.p <- model.summary["length.log", "Pr(>|t|)"]
				male.spring.int <- male.fall.int + model.coefs["factor(season)SPRING"]
				male.spring.int.p <- model.summary["factor(season)SPRING", "Pr(>|t|)"]
				male.spring.slope <- male.fall.slope + model.coefs["length.log:factor(season)SPRING"]
				male.spring.slope.p <- model.summary["length.log:factor(season)SPRING", "Pr(>|t|)"]
				female.fall.int <- male.fall.int + model.coefs["factor(sex)2"]
				female.fall.int.p <- model.summary["factor(sex)2", "Pr(>|t|)"]
				female.fall.slope <- male.fall.slope + model.coefs["length.log:factor(sex)2"]
				female.fall.slope.p <- model.summary["length.log:factor(sex)2", "Pr(>|t|)"]
				female.spring.int <- male.spring.int + model.coefs["factor(sex)2"] + model.coefs["factor(sex)2:factor(season)SPRING"]
				female.spring.int.p <- model.summary["factor(sex)2:factor(season)SPRING",  "Pr(>|t|)"]
				female.spring.slope <- male.spring.slope + model.coefs["length.log:factor(sex)2"] + model.coefs["length.log:factor(sex)2:factor(season)SPRING"]
				female.spring.slope.p <- model.summary["length.log:factor(sex)2:factor(season)SPRING", "Pr(>|t|)"]
				length.log <- rep(log(seq(length.range[1], length.range[2], by=.5)),4)
				sex <- c(rep("1", length(length.log)/2), rep("2", length(length.log)/2))
				season <- rep(c(rep("FALL", length(length.log)/4), rep("SPRING", length(length.log)/4)),2)
				predict.length <- data.frame(length.log, sex, season)
				check.diffs <- plot.lw(this.data, this.model, c("sex", "season"), sp.name, predict.length)
				lw.output <- rbind(lw.output, c(sp.name, sp, ds, 1, "FALL", male.fall.slope, male.fall.slope.p, male.fall.int, male.fall.int.p, male.fall.range[1], male.fall.range[2], check.diffs[1]))
				lw.output <- rbind(lw.output, c(sp.name, sp, ds, 1, "SPRING", male.spring.slope, male.spring.slope.p, male.spring.int, male.spring.int.p, male.spring.range[1], male.spring.range[2],check.diffs[2]))
				lw.output <- rbind(lw.output, c(sp.name, sp, ds, 2, "FALL", female.fall.slope, female.fall.slope.p, female.fall.int, female.fall.int.p, female.fall.range[1], female.fall.range[2],check.diffs[3]))
				lw.output <- rbind(lw.output, c(sp.name, sp, ds, 2, "SPRING", female.spring.slope, female.spring.slope.p, female.spring.int, female.spring.int.p, female.spring.range[1], female.spring.range[2],check.diffs[4]))
			}
			else{
				print(paste("Model parameters not significantly different for", model.definition))
			}
		}
	}
lw.output <- lw.output[!is.na(lw.output$species.code),]
lw.output
}
#Relative Condition:
proc import datafile = "lw_parameters.csv"
 out = LWparams
 dbms = csv
 replace;
 getnames = yes;
run;

data LWparams; set LWparams;
 if LW_SVSPP = 13 then svspp = '013';
 if LW_SVSPP = 15 then svspp = '015';
 if LW_SVSPP = 23 then svspp = '023';
 if LW_SVSPP = 26 then svspp = '026';
 if LW_SVSPP = 28 then svspp = '028';
 if LW_SVSPP = 32 then svspp = '032';
 if LW_SVSPP = 72 then svspp = '072';
 if LW_SVSPP = 73 then svspp = '073';
 if LW_SVSPP = 74 then svspp = '074';
 if LW_SVSPP = 75 then svspp = '075';
 if LW_SVSPP = 76 then svspp = '076';
 if LW_SVSPP = 77 then svspp = '077';
 if LW_SVSPP = 78 then svspp = '078';
 if LW_SVSPP ge 100 then svspp = LW_SVSPP;
 if sexMF = 'M' then sex = '1';
 if sexMF = 'F' then sex = '2';
 if sexMF = ' ' then sex = '0';
 if EXPONENT_FALL = . then EXPONENT_FALL=SEASONLESS_EXPONENT;
 if EXPONENT_SPRING = . then EXPONENT_SPRING=SEASONLESS_EXPONENT;
 if COEFFICIENT_FALL = . then COEFFICIENT_FALL=SEASONLESS_COEFFICIENT;
 if COEFFICIENT_SPRING = . then COEFFICIENT_SPRING=SEASONLESS_COEFFICIENT;

proc sort data=LWparams;
 by svspp sex;

proc sort data=lenwt;
 by svspp sex;

data lwdatpar (keep =cruise6 stratum tow station year month day time beglat beglon setdepth
  surftemp bottemp svspp sex length age maturity indid indwt stom_volume stom_wgt expcatchwt expcatchnum
  COEFFICIENT_SPRING EXPONENT_SPRING COEFFICIENT_FALL EXPONENT_FALL SEASONLESS_COEFFICIENT 
  SEASONLESS_EXPONENT);
 merge lenwt (in=d) LWparams (in=p);
 by svspp sex;

data sortlw; set lwdatpar;
 proc sort; by svspp sex year;

data lwdata; set sortlw;
 if indwt = . then delete;
 if length = . then delete;
 if indwt >0;
 svspp1 = svspp*1;
 indwtg=indwt*1000.0;
 cond=indwtg/(length**3);
 if EXPONENT_FALL gt 0 then predwt = (exp(COEFFICIENT_FALL))*length**EXPONENT_FALL;
 if EXPONENT_FALL = . then predwt = (exp(SEASONLESS_COEFFICIENT))*length**SEASONLESS_EXPONENT;
  if EXPONENT_FALL gt 0 then predwtPK = exp(COEFFICIENT_FALL+(EXPONENT_FALL*log(length)));
 if EXPONENT_FALL = . then predwtPK = exp(SEASONLESS_COEFFICIENT+(SEASONLESS_EXPONENT*log(length)));

***Relative condition;
 RelWt = indwt/predwt*100;

proc sort data=lwdata;
 by svspp1 sex year;
run;

```

### Further metadata
#### Public availability statement
The NEFSC fall bottom trawl survey data is available through InPort:
https://inport.nmfs.noaa.gov/inport/item/22560


#### Point on contact - Whom should be contacted?
Laurel Smith

#### Data steward (can be POC; in charge of keeping data updated)
Laurel Smith


### References