wolf.ageK.PVA2

			JFM  5/4/2020



Template for conducting PVA for WDFW wolf plan (WDFW 2011)

Simulations include wolf age structure, approximated from NW MT values.

Include options for population ceiling and floor.



> # WA wolf plan PVA

> # demographic data from NWMT, appendix G

> # WA Wolf Plan, Appendix G: Stage matrix, NW MT



R version 2.11.1 (2010-05-31)



> n.0 <- c(0, 0, 0, 31)

> n.0

[1]  0  0  0 31

> 

> m.wolf <- matrix(c(0, 0.35, 1.04, 1.04,

+                     0.81, 0, 0, 0,

+                     0, 0.52, 0, 0,

+                     0, 0.2, 0.72, 0.72),

+                  byrow=T, nrow=4, ncol=4)

> m.wolf

     [,1] [,2] [,3] [,4]

[1,] 0.00 0.35 1.04 1.04

[2,] 0.81 0.00 0.00 0.00

[3,] 0.00 0.52 0.00 0.00

[4,] 0.00 0.20 0.72 0.72

> n.1 <- m.wolf %*% n.0

> n.1

      [,1]

[1,] 32.24

[2,]  0.00

[3,]  0.00

[4,] 22.32

> n.2 <- m.wolf %*% n.1

> n.2

        [,1]

[1,] 23.2128

[2,] 26.1144

[3,]  0.0000

[4,] 16.0704

> e.wolf <- eigen(m.wolf)

> e.wolf$values[1]

[1] 1.221678+0i

> e.wolf$vectors[,1]

[1] 0.7129652+0i 0.4727119+0i 0.2012070+0i 0.4772209+0i

> # Determine fraction in each age class @ stable age distribution

> e.wolf$vectors[,1]

[1] 0.7129652+0i 0.4727119+0i 0.2012070+0i 0.4772209+0i

> sum(e.wolf$vectors[,1])

[1] 1.864105+0i

> e.wolf$vectors[,1] / sum(e.wolf$vectors[,1])

[1] 0.3824705+0i 0.2535865+0i 0.1079376+0i 0.2560054+0i



> # 38% are pups!



> # Std.Deviations for parameter distributions, as described in project description;

> #    from Maletzke and Wielgus(2011) and derived from Carroll et al.(2006) 

> sd.wolf <- c(0.13, 0.38, 0.16, 0.16, 0.12, 0.04)

> sd.wolf

[1] 0.13 0.38 0.16 0.16 0.12 0.04



> # Start with 50 adult wolves = 25 adult females wolves

> w.0 <- c(0, 0, 0, 25)

> # Test simulation function: single simulation for 10 years; 

> #    ceiling=228 females (=456 total), floor=5 (10 total)

> age.test <- age.sim.var.normK(m.wolf, sd.wolf, w.0, 228, 5, 10)

> length(age.test[1,])

[1] 4

> age.test

          [,1]     [,2]     [,3]     [,4]

 [1,] 38.00000  0.00000  0.00000 22.00000

 [2,] 34.00000 32.00000  0.00000 19.00000

 [3,] 40.00000 29.00000 19.00000 28.00000

 [4,] 75.00000 34.00000 17.00000 50.00000

 [5,] 94.85496 55.69466 17.40458 60.04580

 [6,] 90.65060 55.62651 22.66265 59.06024

 [7,] 92.41420 52.61538 22.26036 60.71006

 [8,] 93.08876 53.28994 20.91124 60.71006

 [9,] 92.41420 53.96450 21.58580 60.03550

[10,] 92.68843 53.44807 21.64985 60.21365

> 

> # Run single simulation for 100 years

> age.test2 <- age.sim.var.normK(m.wolf, sd.wolf, w.0, 228, 5, 100)

> length(age.test2[1,])

[1] 4

> # Show population at end of simulation

> sum(age.test2[100,])

[1] 228

> # Simulation reached population ceiling

> # Show last 10 years of simulation

> age.test2[90:100,]

          [,1]     [,2]     [,3]     [,4]

 [1,] 78.32653 65.14286 31.02041 53.51020

 [2,] 78.08219 65.58904 31.23288 53.09589

 [3,] 78.08219 64.80822 32.01370 53.09589

 [4,] 78.08219 64.80822 31.23288 53.87671

 [5,] 78.59386 64.58703 31.12628 53.69283

 [6,] 78.32653 65.14286 31.02041 53.51020

 [7,] 78.08219 65.58904 31.23288 53.09589

 [8,] 78.08219 64.80822 32.01370 53.09589

 [9,] 78.08219 64.80822 31.23288 53.87671

[10,] 78.59386 64.58703 31.12628 53.69283

[11,] 78.32653 65.14286 31.02041 53.51020

> 

> yr <- c(1:100)

> plot(yr,age.test2[,4])



> # Test PVA function w/ only 10 simulations, before scaling up.

> pva.ageK.test <- sim.age.pva.normK(m.wolf, sd.wolf, w.0, 228, 5, 20, 10)

> length(pva.ageK.test)

[1] 44

> pva.ageK.test

            [,1]     [,2]     [,3]     [,4]

              NA       NA       NA       NA

sim.end 91.89354 58.95057 24.27376 52.88213

sim.end 69.00000 47.00000 32.00000 63.00000

sim.end 83.33793 59.75172 24.37241 60.53793

sim.end 89.52147 62.94479 26.57669 48.95706

sim.end 23.00000 18.00000  3.00000 16.00000

sim.end  1.00000  1.00000  0.00000  0.00000

sim.end 12.00000 10.00000  4.00000  8.00000

sim.end 98.71264 52.41379 27.08046 49.79310

sim.end 12.00000  7.00000  2.00000 12.00000

sim.end 84.07973 49.99336 23.48173 70.44518

> 

> # Start with 10 adult females, use floor of only 5, ceiling=228 females

> w.1 <- c(0, 0, 0, 10)

> pva.ageK.test3 <- sim.age.pva.normK(m.wolf, sd.wolf, w.1, 228, 5, 100, 1000)

> # Show results for similations 10-20:

> pva.ageK.test3[10:20,]

            [,1]     [,2]     [,3]     [,4]

sim.end  1.00000  1.00000  0.00000  1.00000

sim.end 97.47826 60.30435 23.13043 47.08696

sim.end  1.00000  0.00000  0.00000  1.00000

sim.end 72.83333 60.95833 27.70833 66.50000

sim.end 89.45833 67.29167 19.00000 52.25000

sim.end  1.00000  1.00000  0.00000  1.00000

sim.end 90.51685 58.06742 25.61798 53.79775

sim.end 94.02062 47.79381 26.63918 59.54639

sim.end  1.00000  1.00000  0.00000  2.00000

sim.end  5.00000  3.00000  1.00000  5.00000

sim.end 89.27324 50.09577 31.47042 57.16056

> 

> # Calculate population totals @ end of each simulation

> w.sum3 <- rep(0,1000)

> for (i in c(1:1000)) {w.sum3[i] <- sum(pva.ageK.test3[i,])}

> hist(w.sum3, xlim=c(0,250))

> 

> pva.ageK.test4 <- sim.age.pva.normK(m.wolf, sd.wolf, w.1, 228, 5, 20, 1000)

> # Calculate population totals @ end of each simulation

> w.sum4 <- rep(0,1000)

> for (i in c(1:1000)) {w.sum4[i] <- sum(pva.ageK.test4[i,])}

> hist(w.sum4, xlim=c(0,250), xlab="Total abundance in year 100", ylab="Number of simulations")

> 







# ~~~~~~~~~~~~~~~~~~ R commands only ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~



m.wolf <- matrix(c(0, 0.35, 1.04, 1.04,

                   0.81, 0, 0, 0,

                   0, 0.52, 0, 0,

                   0, 0.2, 0.72, 0.72),

                byrow=T, nrow=4, ncol=4)

m.wolf



# Std.Deviations for parameter distributions, as described in project description;

#    from Maletzke and Wielgus(2011) and derived from Carroll et al.(2006) 

sd.wolf <- c(0.13, 0.38, 0.16, 0.16, 0.12, 0.04)

sd.wolf



# Start with 50 adult wolves = 25 adult females wolves

w.0 <- c(0, 0, 0, 25)

# Test simulation function: single simulation for 10 years; 

#    ceiling=228 females (=456 total), floor=5 (10 total)

age.test <- age.sim.var.normK(m.wolf, sd.wolf, w.0, 228, 5, 10)

length(age.test[1,])

age.test



# Run single simulation for 100 years

age.test2 <- age.sim.var.normK(m.wolf, sd.wolf, w.0, 228, 5, 100)

length(age.test2[1,])

# Show population at end of simulation

sum(age.test2[100,])

# Simulation reached population ceiling

# Show last 10 years of simulation

age.test2[90:100,]



yr <- c(1:100)

plot(yr,age.test2[,4])



# Test PVA function w/ only 10 simulations, before scaling up.

pva.ageK.test <- sim.age.pva.normK(m.wolf, sd.wolf, w.0, 228, 5, 20, 10)

length(pva.ageK.test)

pva.ageK.test



# Start with 10 adult females, use floor of only 5, ceiling=228 females

w.1 <- c(0, 0, 0, 10)

pva.ageK.test3 <- sim.age.pva.normK(m.wolf, sd.wolf, w.1, 228, 5, 100, 1000)

# Show results for similations 10-20:

pva.ageK.test3[10:20,]



# Calculate population totals @ end of each simulation

w.sum3 <- rep(0,1000)

for (i in c(1:1000)) {w.sum3[i] <- sum(pva.ageK.test3[i,])}

hist(w.sum3, xlim=c(0,250))



pva.ageK.test4 <- sim.age.pva.normK(m.wolf, sd.wolf, w.1, 228, 5, 20, 1000)

# Calculate population totals @ end of each simulation

w.sum4 <- rep(0,1000)

for (i in c(1:1000)) {w.sum4[i] <- sum(pva.ageK.test4[i,])}

hist(w.sum4, xlim=c(0,250), xlab="Total abundance in year 100", ylab="Number of simulations")



# ~~~~~~~~~~~~~~~~~~ PVA functions; need both ~~~~~~~~~~~~~~~~~~~~~~~



age.sim.var.normK <- function(m1, mdev, nstart, Klim, floor, it.len) {

    # simulates trajectory of age or stage-structured population,

    #   where fecundities of 3rd and 4th ages/stages normally distributed

    #   and survival rates normally distributed

    # Function arguments:

    #	m1: Age or stage-structured population matrix

    #   mdev: standard deviations: [F2, F3&4, P134, P23, P24] 

    #	nstart:   Starting population, vector of age class abundances

    #   Klim: population ceiling, summed over all age classes

    #   floor: population floor; when sum(nt)< floor, simulation ends

    # 	it.len:   length of iteration (e.g., 20, for 20 yr)

      zero <- 0.0000000000000000

      m1.mean <- m1

      m1.var <- m1.mean

      m1.var[1,2] <- rnorm(1, m1[1,2], mdev[1])  # age2 fecundity

      m1.var[1,3] <- rnorm(1, m1[1,3], mdev[2])  # age3 fecundity 

      m1.var[1,4] <- rnorm(1, m1[1,4], mdev[2])  # age4 fecundity

      m1.var[2,1] <- rnorm(1, m1[2,1], mdev[3])  # age1 survival 

      m1.var[4,3] <- rnorm(1, m1[4,3], mdev[3])  # age3 survival

      m1.var[4,4] <- rnorm(1, m1[4,4], mdev[3])  # age4 survival

      m1.var[3,2] <- rnorm(1, m1[3,2], mdev[4])  # survival juv->disp

      m1.var[4,2] <- rnorm(1, m1[4,2], mdev[5])  # survival juv->adult

      nt <- nstart 

      nt1 <- nstart

      ntraj <- nstart

      dim <- length(nstart)

      ntraj <- matrix(nrow=it.len, ncol=dim)

      j <- 0	

      while (j < it.len) {

          j <- j + 1

      nt <- m1.var %*% nt1

      nt <- trunc(nt)

      if (sum(nt) > Klim) {

          nt <- nt * Klim/sum(nt)

          }

      if (sum(nt) < floor) { 

    #                           nt <- zero

          j <- it.len

          }

      nt1 <- nt

      ntraj[j,] <- nt

      }

      ntraj

}



sim.age.pva.normK <- function(m1, mdev, nstart, Klim, floor, it.len, sim.num) { 

    # Performs "population viability analysis" for age-structured population.

    #     Created for population with 4 age/stage classes; easily modified 

    #        to other matrix dimensions. 

    #     Fecundities and survival rates normally distributed.  

    # Output: An array of the population age distributions at the end of each

    #     simulation.

    # Each row i contains number in each age class at end of the ith

    #     simulation.

    # A row containing all zeros implies extinction for that simulation.

    # More sophisticated versions would treat all matrix elements as variable.

    #     That project is left as an exercise. 

    # Function arguments: 

    #	  m1: Matrix of age/stage-specific fecundities and survival rates, 

    #         must be square matrix.

    #	  nstart:   starting population age distribution, 

    #               must be same dimension as matrix. 

    #     it.len:   length of iteration (e.g., 20, for 20 yr)

    #     sim.num:  number of simulations to run

    # Calls function "age.sim.var.normK", which simulates pop growth age 

    #   structure with normally distributed Fx and Px 

        col <- length(nstart)

 	sim.i <- NULL

	sim.res <- matrix(ncol=col)

	i <- 0

        while (i < sim.num) {

	    i <- i + 1

	    sim.end <- nstart

            sim.i <- age.sim.var.normK(m1, mdev, nstart, Klim, floor, it.len)

            sim.end <- sim.i[it.len,]

	    sim.res <- rbind(sim.res, sim.end)

        }

        sim.res

}