w <- pp / sum(pp)
denom <- log(pp)
}
ki2[j] <-  0.25 +  sum( w * (log(TE[j])/denom) )
}
B <- (mean.M)^ki2                                         # biomass proportional to body size
k <- summary(lm( c(log(B)+0.001) ~ log(mean.M)))$coef[2,1]
harvest.B <- B
harvest.B[fished_sp] <- (1-df$remove.bm) * B[fished_sp]
harvest.k <- summary(lm(log(harvest.B) ~ log(harvest.mean.M)))$coef[2,1]
# harvest.TE <- sTE * harvest.mean.M^(-0.07)                 # transfer efficiency dependent on body size (Woodson et al 2018)
# harvest.PPMR <- exp(log(harvest.mean.M) / (TL-1+E) )       # transfer efficiency dependent on body size (Woodson et al 2018)
# harvest.ki <-  0.25 + log(harvest.TE) / log(harvest.PPMR)  # scaling coefficient
df[i, 1:9] <- c( sp=sp,
sTE=sTE[i],
meanTE=mean(TE),
k=k,
harvest.k=harvest.k,
PPMR=mean(PPMR),
maxPPMR=max(PPMR),
medPPMR=median(PPMR),
maxTL=max(TL) )
}
}
i
df
db <- list()
for(i in 1:n){
C <- sp^(0.35-1)                         # connectance
FM <- trophic::niche(sp, C)              # feeding matrix
TL <- TrophInd(FM)[,1]                   # trophic levels
E <- rnorm(sp, mean=0, sd=0.1)          # noise to BS
PPMR <- exp(1) * TL^4.8                 # original
# PPMR <- 10^3 * TL^rnorm(sp, 0, 0.1)
#PPMR <- 10^3
# plot(PPMR, TL)
M0 <- 1
mean.M <- M0*PPMR^(TL-1+E)
sd.M <- 0.1 * mean.M                       # sd of body sizes increases with the mean
q <- runif(sp, 0.6, 1)                                   # exponent of search volume (Andersen et al 2009)
TE <- sTE[i] * TL^(0.5-q)                   # transfer efficiency (Andersen et al 2009)
#  plot(TL, TE)
# size-selective harvest
probs.M <- mean.M/sum(mean.M)              # prob of being fished based on size
fished_sp <- which(mean.M %in% sample(mean.M, floor(sp*frac), prob=probs.M))  # sample larger species
remove <- runif(sp, 0.3, 1)
harvest <- lapply(fished_sp, function(x) {
dist <-  rnorm( 5000, mean = mean.M[[x]], sd = sd.M[[x]] )
probs <- pnorm(dist, mean = mean.M[[x]], sd = sd.M[[x]])
harvest.dist <- dist[-which(dist %in% sample(dist, floor(remove[[x]]*5000), prob=probs))]
new.mean.M <- mean(harvest.dist)
sd.new.M <- sd(harvest.dist)
remove.bm <- (sum(dist)  - sum(harvest.dist)) / sum(dist)
return(data.frame(sp=x, mean.M=mean.M[[x]], new.mean.M=new.mean.M, remove.bm=remove.bm, remove.ab=remove[[x]]))
})
df0 <- suppressMessages({ melt(harvest) })
df <- tidyr::spread(df0, variable, value)
harvest.mean.M <-  mean.M
harvest.mean.M[fished_sp] <- df$new.mean.M
if( any(TE>1) ) {
return(NA)
}else{
PPMRs <- FM * matrix(rep(mean.M, sp), nrow=sp, byrow=T) / matrix(rep(mean.M, sp), nrow=sp, byrow=F)
ki <-  0.25 + log(TE) / log(PPMR)  # scaling coefficient
ki2 <- NULL
for(j in 1:sp){
pp <- subset(PPMRs[,j], PPMRs[,j]>0 )
pp[which(pp==1)] <- pp[which(pp==1)]*0.9
if(length(pp)==0){
w <- 1
denom <- log(PPMR[j] )
} else {
w <- pp / sum(pp)
denom <- log(pp)
}
ki2[j] <-  0.25 +  sum( w * (log(TE[j])/denom) )
}
B <- (mean.M)^ki2                                         # biomass proportional to body size
k <- summary(lm( c(log(B)+0.001) ~ log(mean.M)))$coef[2,1]
harvest.B <- B
harvest.B[fished_sp] <- (1-df$remove.bm) * B[fished_sp]
harvest.k <- summary(lm(log(harvest.B) ~ log(harvest.mean.M)))$coef[2,1]
# harvest.TE <- sTE * harvest.mean.M^(-0.07)                 # transfer efficiency dependent on body size (Woodson et al 2018)
# harvest.PPMR <- exp(log(harvest.mean.M) / (TL-1+E) )       # transfer efficiency dependent on body size (Woodson et al 2018)
# harvest.ki <-  0.25 + log(harvest.TE) / log(harvest.PPMR)  # scaling coefficient
db[i, 1:9] <- c( sp=sp,
sTE=sTE[i],
meanTE=mean(TE),
k=k,
harvest.k=harvest.k,
PPMR=mean(PPMR),
maxPPMR=max(PPMR),
medPPMR=median(PPMR),
maxTL=max(TL) )
}
}
db[i,] <- c( sp=sp,
sTE=sTE[i],
meanTE=mean(TE),
k=k,
harvest.k=harvest.k,
PPMR=mean(PPMR),
maxPPMR=max(PPMR),
medPPMR=median(PPMR),
maxTL=max(TL) )
db[[i]] <- c( sp=sp,
sTE=sTE[i],
meanTE=mean(TE),
k=k,
harvest.k=harvest.k,
PPMR=mean(PPMR),
maxPPMR=max(PPMR),
medPPMR=median(PPMR),
maxTL=max(TL) )
db <- list()
for(i in 1:n){
C <- sp^(0.35-1)                         # connectance
FM <- trophic::niche(sp, C)              # feeding matrix
TL <- TrophInd(FM)[,1]                   # trophic levels
E <- rnorm(sp, mean=0, sd=0.1)          # noise to BS
PPMR <- exp(1) * TL^4.8                 # original
# PPMR <- 10^3 * TL^rnorm(sp, 0, 0.1)
#PPMR <- 10^3
# plot(PPMR, TL)
M0 <- 1
mean.M <- M0*PPMR^(TL-1+E)
sd.M <- 0.1 * mean.M                       # sd of body sizes increases with the mean
q <- runif(sp, 0.6, 1)                                   # exponent of search volume (Andersen et al 2009)
TE <- sTE[i] * TL^(0.5-q)                   # transfer efficiency (Andersen et al 2009)
#  plot(TL, TE)
# size-selective harvest
probs.M <- mean.M/sum(mean.M)              # prob of being fished based on size
fished_sp <- which(mean.M %in% sample(mean.M, floor(sp*frac), prob=probs.M))  # sample larger species
remove <- runif(sp, 0.3, 1)
harvest <- lapply(fished_sp, function(x) {
dist <-  rnorm( 5000, mean = mean.M[[x]], sd = sd.M[[x]] )
probs <- pnorm(dist, mean = mean.M[[x]], sd = sd.M[[x]])
harvest.dist <- dist[-which(dist %in% sample(dist, floor(remove[[x]]*5000), prob=probs))]
new.mean.M <- mean(harvest.dist)
sd.new.M <- sd(harvest.dist)
remove.bm <- (sum(dist)  - sum(harvest.dist)) / sum(dist)
return(data.frame(sp=x, mean.M=mean.M[[x]], new.mean.M=new.mean.M, remove.bm=remove.bm, remove.ab=remove[[x]]))
})
df0 <- suppressMessages({ melt(harvest) })
df <- tidyr::spread(df0, variable, value)
harvest.mean.M <-  mean.M
harvest.mean.M[fished_sp] <- df$new.mean.M
if( any(TE>1) ) {
return(NA)
}else{
PPMRs <- FM * matrix(rep(mean.M, sp), nrow=sp, byrow=T) / matrix(rep(mean.M, sp), nrow=sp, byrow=F)
ki <-  0.25 + log(TE) / log(PPMR)  # scaling coefficient
ki2 <- NULL
for(j in 1:sp){
pp <- subset(PPMRs[,j], PPMRs[,j]>0 )
pp[which(pp==1)] <- pp[which(pp==1)]*0.9
if(length(pp)==0){
w <- 1
denom <- log(PPMR[j] )
} else {
w <- pp / sum(pp)
denom <- log(pp)
}
ki2[j] <-  0.25 +  sum( w * (log(TE[j])/denom) )
}
B <- (mean.M)^ki2                                         # biomass proportional to body size
k <- summary(lm( c(log(B)+0.001) ~ log(mean.M)))$coef[2,1]
harvest.B <- B
harvest.B[fished_sp] <- (1-df$remove.bm) * B[fished_sp]
harvest.k <- summary(lm(log(harvest.B) ~ log(harvest.mean.M)))$coef[2,1]
# harvest.TE <- sTE * harvest.mean.M^(-0.07)                 # transfer efficiency dependent on body size (Woodson et al 2018)
# harvest.PPMR <- exp(log(harvest.mean.M) / (TL-1+E) )       # transfer efficiency dependent on body size (Woodson et al 2018)
# harvest.ki <-  0.25 + log(harvest.TE) / log(harvest.PPMR)  # scaling coefficient
db[[i]] <- c( sp=sp,
sTE=sTE[i],
meanTE=mean(TE),
k=k,
harvest.k=harvest.k,
PPMR=mean(PPMR),
maxPPMR=max(PPMR),
medPPMR=median(PPMR),
maxTL=max(TL) )
}
}
i
db
n
db <- list()
for(i in 1:n){
C <- sp^(0.35-1)                         # connectance
FM <- trophic::niche(sp, C)              # feeding matrix
TL <- TrophInd(FM)[,1]                   # trophic levels
E <- rnorm(sp, mean=0, sd=0.1)          # noise to BS
PPMR <- exp(1) * TL^4.8                 # original
# PPMR <- 10^3 * TL^rnorm(sp, 0, 0.1)
#PPMR <- 10^3
# plot(PPMR, TL)
M0 <- 1
mean.M <- M0*PPMR^(TL-1+E)
sd.M <- 0.1 * mean.M                       # sd of body sizes increases with the mean
q <- runif(sp, 0.6, 1)                                   # exponent of search volume (Andersen et al 2009)
TE <- sTE[i] * TL^(0.5-q)                   # transfer efficiency (Andersen et al 2009)
#  plot(TL, TE)
# size-selective harvest
probs.M <- mean.M/sum(mean.M)              # prob of being fished based on size
fished_sp <- which(mean.M %in% sample(mean.M, floor(sp*frac), prob=probs.M))  # sample larger species
remove <- runif(sp, 0.3, 1)
harvest <- lapply(fished_sp, function(x) {
dist <-  rnorm( 5000, mean = mean.M[[x]], sd = sd.M[[x]] )
probs <- pnorm(dist, mean = mean.M[[x]], sd = sd.M[[x]])
harvest.dist <- dist[-which(dist %in% sample(dist, floor(remove[[x]]*5000), prob=probs))]
new.mean.M <- mean(harvest.dist)
sd.new.M <- sd(harvest.dist)
remove.bm <- (sum(dist)  - sum(harvest.dist)) / sum(dist)
return(data.frame(sp=x, mean.M=mean.M[[x]], new.mean.M=new.mean.M, remove.bm=remove.bm, remove.ab=remove[[x]]))
})
df0 <- suppressMessages({ melt(harvest) })
df <- tidyr::spread(df0, variable, value)
harvest.mean.M <-  mean.M
harvest.mean.M[fished_sp] <- df$new.mean.M
if( any(TE>1) ) {
return(NA)
}else{
PPMRs <- FM * matrix(rep(mean.M, sp), nrow=sp, byrow=T) / matrix(rep(mean.M, sp), nrow=sp, byrow=F)
ki <-  0.25 + log(TE) / log(PPMR)  # scaling coefficient
ki2 <- NULL
for(j in 1:sp){
pp <- subset(PPMRs[,j], PPMRs[,j]>0 )
pp[which(pp==1)] <- pp[which(pp==1)]*0.9
if(length(pp)==0){
w <- 1
denom <- log(PPMR[j] )
} else {
w <- pp / sum(pp)
denom <- log(pp)
}
ki2[j] <-  0.25 +  sum( w * (log(TE[j])/denom) )
}
B <- (mean.M)^ki2                                         # biomass proportional to body size
k <- summary(lm( c(log(B)+0.001) ~ log(mean.M)))$coef[2,1]
harvest.B <- B
harvest.B[fished_sp] <- (1-df$remove.bm) * B[fished_sp]
harvest.k <- summary(lm(log(harvest.B) ~ log(harvest.mean.M)))$coef[2,1]
# harvest.TE <- sTE * harvest.mean.M^(-0.07)                 # transfer efficiency dependent on body size (Woodson et al 2018)
# harvest.PPMR <- exp(log(harvest.mean.M) / (TL-1+E) )       # transfer efficiency dependent on body size (Woodson et al 2018)
# harvest.ki <-  0.25 + log(harvest.TE) / log(harvest.PPMR)  # scaling coefficient
db[[i]] <- c( sp=sp,
sTE=sTE[i],
meanTE=mean(TE),
k=k,
harvest.k=harvest.k,
PPMR=mean(PPMR),
maxPPMR=max(PPMR),
medPPMR=median(PPMR),
maxTL=max(TL) )
}
}
db
C <- sp^(0.35-1)                         # connectance
FM <- trophic::niche(sp, C)              # feeding matrix
TL <- TrophInd(FM)[,1]                   # trophic levels
E <- rnorm(sp, mean=0, sd=0.1)          # noise to BS
PPMR <- exp(1) * TL^4.8                 # original
M0 <- 1
mean.M <- M0*PPMR^(TL-1+E)
sd.M <- 0.1 * mean.M                       # sd of body sizes increases with the mean
q <- runif(sp, 0.6, 1)                                   # exponent of search volume (Andersen et al 2009)
TE <- sTE[i] * TL^(0.5-q)                   # transfer efficiency (Andersen et al 2009)
# size-selective harvest
probs.M <- mean.M/sum(mean.M)              # prob of being fished based on size
fished_sp <- which(mean.M %in% sample(mean.M, floor(sp*frac), prob=probs.M))  # sample larger species
remove <- runif(sp, 0.3, 1)
harvest <- lapply(fished_sp, function(x) {
dist <-  rnorm( 5000, mean = mean.M[[x]], sd = sd.M[[x]] )
probs <- pnorm(dist, mean = mean.M[[x]], sd = sd.M[[x]])
harvest.dist <- dist[-which(dist %in% sample(dist, floor(remove[[x]]*5000), prob=probs))]
new.mean.M <- mean(harvest.dist)
sd.new.M <- sd(harvest.dist)
remove.bm <- (sum(dist)  - sum(harvest.dist)) / sum(dist)
return(data.frame(sp=x, mean.M=mean.M[[x]], new.mean.M=new.mean.M, remove.bm=remove.bm, remove.ab=remove[[x]]))
})
df0 <- suppressMessages({ melt(harvest) })
df <- tidyr::spread(df0, variable, value)
harvest.mean.M <-  mean.M
harvest.mean.M[fished_sp] <- df$new.mean.M
if( any(TE>1) ) {
return(NA)
}else{
PPMRs <- FM * matrix(rep(mean.M, sp), nrow=sp, byrow=T) / matrix(rep(mean.M, sp), nrow=sp, byrow=F)
ki <-  0.25 + log(TE) / log(PPMR)  # scaling coefficient
ki2 <- NULL
for(j in 1:sp){
pp <- subset(PPMRs[,j], PPMRs[,j]>0 )
pp[which(pp==1)] <- pp[which(pp==1)]*0.9
if(length(pp)==0){
w <- 1
denom <- log(PPMR[j] )
} else {
w <- pp / sum(pp)
denom <- log(pp)
}
ki2[j] <-  0.25 +  sum( w * (log(TE[j])/denom) )
}
B <- (mean.M)^ki2                                         # biomass proportional to body size
k <- summary(lm( c(log(B)+0.001) ~ log(mean.M)))$coef[2,1]
harvest.B <- B
harvest.B[fished_sp] <- (1-df$remove.bm) * B[fished_sp]
harvest.k <- summary(lm(log(harvest.B) ~ log(harvest.mean.M)))$coef[2,1]
# harvest.TE <- sTE * harvest.mean.M^(-0.07)                 # transfer efficiency dependent on body size (Woodson et al 2018)
# harvest.PPMR <- exp(log(harvest.mean.M) / (TL-1+E) )       # transfer efficiency dependent on body size (Woodson et al 2018)
# harvest.ki <-  0.25 + log(harvest.TE) / log(harvest.PPMR)  # scaling coefficient
db[[i]] <- c( sp=sp,
sTE=sTE[i],
meanTE=mean(TE),
k=k,
harvest.k=harvest.k,
PPMR=mean(PPMR),
maxPPMR=max(PPMR),
medPPMR=median(PPMR),
maxTL=max(TL) )
}
PPMRs <- FM * matrix(rep(mean.M, sp), nrow=sp, byrow=T) / matrix(rep(mean.M, sp), nrow=sp, byrow=F)
ki <-  0.25 + log(TE) / log(PPMR)  # scaling coefficient
ki2 <- NULL
for(j in 1:sp){
pp <- subset(PPMRs[,j], PPMRs[,j]>0 )
pp[which(pp==1)] <- pp[which(pp==1)]*0.9
if(length(pp)==0){
w <- 1
denom <- log(PPMR[j] )
} else {
w <- pp / sum(pp)
denom <- log(pp)
}
ki2[j] <-  0.25 +  sum( w * (log(TE[j])/denom) )
}
B <- (mean.M)^ki2                                         # biomass proportional to body size
k <- summary(lm( c(log(B)+0.001) ~ log(mean.M)))$coef[2,1]
harvest.B <- B
harvest.B[fished_sp] <- (1-df$remove.bm) * B[fished_sp]
harvest.k <- summary(lm(log(harvest.B) ~ log(harvest.mean.M)))$coef[2,1]
db[[i]] <- c( sp=sp,
sTE=sTE[i],
meanTE=mean(TE),
k=k,
harvest.k=harvest.k,
PPMR=mean(PPMR),
maxPPMR=max(PPMR),
medPPMR=median(PPMR),
maxTL=max(TL) )
db
db[[i]] <- c(      n = i,
sp=sp,
sTE=sTE[i],
meanTE=mean(TE),
k=k,
harvest.k=harvest.k,
PPMR=mean(PPMR),
maxPPMR=max(PPMR),
medPPMR=median(PPMR),
maxTL=max(TL) )
db
db <- list()
for(i in 1:n){
C <- sp^(0.35-1)                         # connectance
FM <- trophic::niche(sp, C)              # feeding matrix
TL <- TrophInd(FM)[,1]                   # trophic levels
E <- rnorm(sp, mean=0, sd=0.1)          # noise to BS
PPMR <- exp(1) * TL^4.8                 # original
# PPMR <- 10^3 * TL^rnorm(sp, 0, 0.1)
#PPMR <- 10^3
# plot(PPMR, TL)
M0 <- 1
mean.M <- M0*PPMR^(TL-1+E)
sd.M <- 0.1 * mean.M                       # sd of body sizes increases with the mean
q <- runif(sp, 0.6, 1)                                   # exponent of search volume (Andersen et al 2009)
TE <- sTE[i] * TL^(0.5-q)                   # transfer efficiency (Andersen et al 2009)
#  plot(TL, TE)
# size-selective harvest
probs.M <- mean.M/sum(mean.M)              # prob of being fished based on size
fished_sp <- which(mean.M %in% sample(mean.M, floor(sp*frac), prob=probs.M))  # sample larger species
remove <- runif(sp, 0.3, 1)
harvest <- lapply(fished_sp, function(x) {
dist <-  rnorm( 5000, mean = mean.M[[x]], sd = sd.M[[x]] )
probs <- pnorm(dist, mean = mean.M[[x]], sd = sd.M[[x]])
harvest.dist <- dist[-which(dist %in% sample(dist, floor(remove[[x]]*5000), prob=probs))]
new.mean.M <- mean(harvest.dist)
sd.new.M <- sd(harvest.dist)
remove.bm <- (sum(dist)  - sum(harvest.dist)) / sum(dist)
return(data.frame(sp=x, mean.M=mean.M[[x]], new.mean.M=new.mean.M, remove.bm=remove.bm, remove.ab=remove[[x]]))
})
df0 <- suppressMessages({ melt(harvest) })
df <- tidyr::spread(df0, variable, value)
harvest.mean.M <-  mean.M
harvest.mean.M[fished_sp] <- df$new.mean.M
if( any(TE>1) ) {
return(NA)
}else{
PPMRs <- FM * matrix(rep(mean.M, sp), nrow=sp, byrow=T) / matrix(rep(mean.M, sp), nrow=sp, byrow=F)
ki <-  0.25 + log(TE) / log(PPMR)  # scaling coefficient
ki2 <- NULL
for(j in 1:sp){
pp <- subset(PPMRs[,j], PPMRs[,j]>0 )
pp[which(pp==1)] <- pp[which(pp==1)]*0.9
if(length(pp)==0){
w <- 1
denom <- log(PPMR[j] )
} else {
w <- pp / sum(pp)
denom <- log(pp)
}
ki2[j] <-  0.25 +  sum( w * (log(TE[j])/denom) )
}
B <- (mean.M)^ki2                                         # biomass proportional to body size
k <- summary(lm( c(log(B)+0.001) ~ log(mean.M)))$coef[2,1]
harvest.B <- B
harvest.B[fished_sp] <- (1-df$remove.bm) * B[fished_sp]
harvest.k <- summary(lm(log(harvest.B) ~ log(harvest.mean.M)))$coef[2,1]
# harvest.TE <- sTE * harvest.mean.M^(-0.07)                 # transfer efficiency dependent on body size (Woodson et al 2018)
# harvest.PPMR <- exp(log(harvest.mean.M) / (TL-1+E) )       # transfer efficiency dependent on body size (Woodson et al 2018)
# harvest.ki <-  0.25 + log(harvest.TE) / log(harvest.PPMR)  # scaling coefficient
db[[i]] <- data.frame(      n = i,
sp=sp,
sTE=sTE[i],
meanTE=mean(TE),
k=k,
harvest.k=harvest.k,
PPMR=mean(PPMR),
maxPPMR=max(PPMR),
medPPMR=median(PPMR),
maxTL=max(TL) )
}
}
db
sTE
sTE <- round(runif(n, 0.05, 1), 1)    # scaling factor of transfer efficiency
sp <- 50                              # number of species
frac <- 0.3                           # fraction of species being fished
n <- 10^3                              # number of simulations
sTE <- round(runif(n, 0.05, 1), 1)    # scaling factor of transfer efficiency
db <- list()
for(i in 1:n){
C <- sp^(0.35-1)                         # connectance
FM <- trophic::niche(sp, C)              # feeding matrix
TL <- TrophInd(FM)[,1]                   # trophic levels
E <- rnorm(sp, mean=0, sd=0.1)          # noise to BS
PPMR <- exp(1) * TL^4.8                 # original
# PPMR <- 10^3 * TL^rnorm(sp, 0, 0.1)
#PPMR <- 10^3
# plot(PPMR, TL)
M0 <- 1
mean.M <- M0*PPMR^(TL-1+E)
sd.M <- 0.1 * mean.M                       # sd of body sizes increases with the mean
q <- runif(sp, 0.6, 1)                                   # exponent of search volume (Andersen et al 2009)
TE <- sTE[i] * TL^(0.5-q)                   # transfer efficiency (Andersen et al 2009)
#  plot(TL, TE)
# size-selective harvest
probs.M <- mean.M/sum(mean.M)              # prob of being fished based on size
fished_sp <- which(mean.M %in% sample(mean.M, floor(sp*frac), prob=probs.M))  # sample larger species
remove <- runif(sp, 0.3, 1)
harvest <- lapply(fished_sp, function(x) {
dist <-  rnorm( 5000, mean = mean.M[[x]], sd = sd.M[[x]] )
probs <- pnorm(dist, mean = mean.M[[x]], sd = sd.M[[x]])
harvest.dist <- dist[-which(dist %in% sample(dist, floor(remove[[x]]*5000), prob=probs))]
new.mean.M <- mean(harvest.dist)
sd.new.M <- sd(harvest.dist)
remove.bm <- (sum(dist)  - sum(harvest.dist)) / sum(dist)
return(data.frame(sp=x, mean.M=mean.M[[x]], new.mean.M=new.mean.M, remove.bm=remove.bm, remove.ab=remove[[x]]))
})
df0 <- suppressMessages({ melt(harvest) })
df <- tidyr::spread(df0, variable, value)
harvest.mean.M <-  mean.M
harvest.mean.M[fished_sp] <- df$new.mean.M
if( any(TE>1) ) {
return(NA)
}else{
PPMRs <- FM * matrix(rep(mean.M, sp), nrow=sp, byrow=T) / matrix(rep(mean.M, sp), nrow=sp, byrow=F)
ki <-  0.25 + log(TE) / log(PPMR)  # scaling coefficient
ki2 <- NULL
for(j in 1:sp){
pp <- subset(PPMRs[,j], PPMRs[,j]>0 )
pp[which(pp==1)] <- pp[which(pp==1)]*0.9
if(length(pp)==0){
w <- 1
denom <- log(PPMR[j] )
} else {
w <- pp / sum(pp)
denom <- log(pp)
}
ki2[j] <-  0.25 +  sum( w * (log(TE[j])/denom) )
}
B <- (mean.M)^ki2                                         # biomass proportional to body size
k <- summary(lm( c(log(B)+0.001) ~ log(mean.M)))$coef[2,1]
harvest.B <- B
harvest.B[fished_sp] <- (1-df$remove.bm) * B[fished_sp]
harvest.k <- summary(lm(log(harvest.B) ~ log(harvest.mean.M)))$coef[2,1]
# harvest.TE <- sTE * harvest.mean.M^(-0.07)                 # transfer efficiency dependent on body size (Woodson et al 2018)
# harvest.PPMR <- exp(log(harvest.mean.M) / (TL-1+E) )       # transfer efficiency dependent on body size (Woodson et al 2018)
# harvest.ki <-  0.25 + log(harvest.TE) / log(harvest.PPMR)  # scaling coefficient
db[[i]] <- data.frame(      n = i,
sp=sp,
sTE=sTE[i],
meanTE=mean(TE),
k=k,
harvest.k=harvest.k,
PPMR=mean(PPMR),
maxPPMR=max(PPMR),
medPPMR=median(PPMR),
maxTL=max(TL) )
}
}
db
