Positive Bernoulli Sequence Model

Density, and random generation for multiple Bernoulli responses where each row in the response matrix has at least one success.

Usage

rposbern(n, nTimePts = 5, pvars = length(xcoeff),
  xcoeff = c(-2, 1, 2), Xmatrix = NULL, cap.effect = 1,
  is.popn = FALSE, link = "logitlink", earg.link = FALSE)
dposbern(x, prob, prob0 = prob, log = FALSE)

Arguments

x

response vector or matrix. Should only have 0 and 1 values, at least two columns, and each row should have at least one 1.

nTimePts

Number of sampling occasions. Called \(\tau\) in posbernoulli.b and posbernoulli.t.

n

number of observations. Usually a single positive integer, else the length of the vector is used. See argument is.popn.

is.popn

Logical. If TRUE then argument n is the population size and what is returned may have substantially less rows than n. That is, if an animal has at least one one in its sequence then it is returned, else that animal is not returned because it never was captured.

Xmatrix

Optional X matrix. If given, the X matrix is not generated internally.

cap.effect

Numeric, the capture effect. Added to the linear predictor if captured previously. A positive or negative value corresponds to a trap-happy and trap-shy effect respectively.

pvars

Number of other numeric covariates that make up the linear predictor. Labelled x1, x2, ..., where the first is an intercept, and the others are independent standard runif random variates. The first pvars elements of xcoeff are used.

xcoeff

The regression coefficients of the linear predictor. These correspond to x1, x2, ..., and the first is for the intercept. The length of xcoeff must be at least pvars.

link, earg.link

The former is used to generate the probabilities for capture at each occasion. Other details at CommonVGAMffArguments.

prob, prob0

Matrix of probabilities for the numerator and denominators respectively. The default does not correspond to the \(M_b\) model since the \(M_b\) model has a denominator which involves the capture history.

log

Logical. Return the logarithm of the answer?

Details

The form of the conditional likelihood is described in posbernoulli.b and/or posbernoulli.t and/or posbernoulli.tb. The denominator is equally shared among the elements of the matrix x.

Value

rposbern returns a data frame with some attributes. The function generates random deviates (\(\tau\) columns labelled y1, y2, ...) for the response. Some indicator columns are also included (those starting with ch are for previous capture history). The default setting corresponds to a \(M_{bh}\) model that has a single trap-happy effect. Covariates x1, x2, ... have the same affect on capture/recapture at every sampling occasion (see the argument parallel.t in, e.g., posbernoulli.tb).

The function dposbern gives the density,

Author

Thomas W. Yee.

Note

The r-type function is experimental only and does not follow the usual conventions of r-type R functions. It may change a lot in the future. The d-type function is more conventional and is less likely to change.

Examples

rposbern(n = 10)
#>    y1 y2 y3 y4 y5 x1         x2        x3 ch1 ch2 ch3 ch4 ch5
#> 1   0  1  1  0  1  1 0.92196792 0.8023485   0   0   1   1   1
#> 2   0  1  0  0  1  1 0.25248029 0.6712682   0   0   1   1   1
#> 3   0  1  1  1  1  1 0.75143276 0.4183216   0   0   1   1   1
#> 4   0  0  1  1  1  1 0.20433855 0.7535079   0   0   0   1   1
#> 5   0  1  0  1  1  1 0.82924332 0.4149254   0   0   1   1   1
#> 6   0  1  0  0  1  1 0.76660860 0.3198663   0   0   1   1   1
#> 7   0  0  1  0  1  1 0.12065501 0.6464812   0   0   0   1   1
#> 8   1  0  1  1  1  1 0.82122233 0.8718459   0   1   1   1   1
#> 9   0  1  1  0  0  1 0.85388981 0.4064510   0   0   1   1   1
#> 10  0  0  0  0  1  1 0.03792608 0.8893126   0   0   0   0   0
attributes(pdata <- rposbern(n = 100))
#> $names
#>  [1] "y1"  "y2"  "y3"  "y4"  "y5"  "x1"  "x2"  "x3"  "ch1" "ch2" "ch3" "ch4"
#> [13] "ch5"
#> 
#> $row.names
#>   [1] "1"   "2"   "3"   "4"   "5"   "6"   "7"   "8"   "9"   "10"  "11"  "12" 
#>  [13] "13"  "14"  "15"  "16"  "17"  "18"  "19"  "20"  "21"  "22"  "23"  "24" 
#>  [25] "25"  "26"  "27"  "28"  "29"  "30"  "31"  "32"  "33"  "34"  "35"  "36" 
#>  [37] "37"  "38"  "39"  "40"  "41"  "42"  "43"  "44"  "45"  "46"  "47"  "48" 
#>  [49] "49"  "50"  "51"  "52"  "53"  "54"  "55"  "56"  "57"  "58"  "59"  "60" 
#>  [61] "61"  "62"  "63"  "64"  "65"  "66"  "67"  "68"  "69"  "70"  "71"  "72" 
#>  [73] "73"  "74"  "75"  "76"  "77"  "78"  "79"  "80"  "81"  "82"  "83"  "84" 
#>  [85] "85"  "86"  "87"  "88"  "89"  "90"  "91"  "92"  "93"  "94"  "95"  "96" 
#>  [97] "97"  "98"  "99"  "100"
#> 
#> $class
#> [1] "data.frame"
#> 
#> $pvars
#> [1] 3
#> 
#> $nTimePts
#> [1] 5
#> 
#> $cap.effect
#> [1] 1
#> 
#> $is.popn
#> [1] FALSE
#> 
#> $n
#> [1] 100
#> 
M.bh <- vglm(cbind(y1, y2, y3, y4, y5) ~ x2 + x3,
             posbernoulli.b(I2 = FALSE), pdata, trace = TRUE)
#> Iteration 1: loglikelihood = -313.79493
#> Iteration 2: loglikelihood = -313.40993
#> Iteration 3: loglikelihood = -313.40766
#> Iteration 4: loglikelihood = -313.40764
#> Iteration 5: loglikelihood = -313.40764
constraints(M.bh)
#> $`(Intercept)`
#>      [,1] [,2]
#> [1,]    0    1
#> [2,]    1    1
#> 
#> $x2
#>      [,1]
#> [1,]    1
#> [2,]    1
#> 
#> $x3
#>      [,1]
#> [1,]    1
#> [2,]    1
#> 
summary(M.bh)
#> 
#> Call:
#> vglm(formula = cbind(y1, y2, y3, y4, y5) ~ x2 + x3, family = posbernoulli.b(I2 = FALSE), 
#>     data = pdata, trace = TRUE)
#> 
#> Coefficients: 
#>               Estimate Std. Error z value Pr(>|z|)    
#> (Intercept):1   0.9164     0.2419   3.787 0.000152 ***
#> (Intercept):2  -1.7556     0.4126  -4.255 2.09e-05 ***
#> x2              1.0972     0.3884   2.825 0.004730 ** 
#> x3              1.6559     0.4138   4.002 6.28e-05 ***
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> Names of linear predictors: logitlink(pcapture), logitlink(precapture)
#> 
#> Log-likelihood: -313.4076 on 196 degrees of freedom
#> 
#> Number of Fisher scoring iterations: 5 
#> 
#> 
#> Estimate of N:  109.527 
#> 
#> Std. Error of N:  5.523 
#> 
#> Approximate 95 percent confidence interval for N:
#> 98.7 120.35