Confidence limits from covariates

[mFidino]

We recently completed a single-species multiple season occupancy analysis for p, psi, and gamma in PRESENCE and used a number of continuous covariates that we thought would influence these parameters. After the analysis we would like to plot out how individual covariates in the top performing model influenced psi, but are having problems doing so as the top model contains multiple covariates (e.g. neighborhood age and housing density). While it is simple enough to use the logistic link function with the beta estimates to create these probabilities the issue we are having is generating confidence limits for each individual covariate to plot along with them. After looking over phidot it either has to do with the variance-covariance matrix / delta method or creating odds-ratios, but we are uncertain of exactly what to do. We would greatly appreciate it if someone could point us in the right direction.

[darryl]

If you have a logistic regression equation of the form:
logit(psi) = a1 + a2*x + a3*y, for covariate values x and y, create a vector of values:
z=[1 x y], then use matrix multiplication to determine Var(logit(psi)) = z V t(z), where V is the portion of the variance-covariance matrix associated with the psi coefficients and t(z) is the transpose of z. You can then calculate the corresponding SE, calculate limits on the logit scale, then back transform the limits onto the real scale. This probably sounds worse than it really is, you can set it up in Excel without a whole lot of difficulty. To get PRESENCE to give you the variance-covariance matrix, you need to select it as an option when you're running the model.

Good luck.
Darryl

[mFidino]

I will give it a go. Thanks!

[mFidino]

So I had some time today to sit down and have a crack at this but I still feel as if I am doing something a little incorrectly. I created the vector z for each covariate value in a larger matrix(e.g. if I have 10 sites the resulting matrix had dimensions 10 x 3).

From here I tried matrix multiplying each row of the larger z matrix into the particular portion of the variance-covariance matrix (which had dimensions of [3 x 3] as there are three beta coefficients related to psi). For example, for the first of the covariate values:

[1 x 3] %*% [3 x 3] returns a matrix with dimensions [3 x 1].

However, here you cannot multiply this matrix by t(z), which is the transpose of the [1 x 3] matrix from the first bit of matrix math because you would be multiplying [3 X 1] %*% [3 x 1], which you cannot do.

I am likely doing something incorrect here, but am having a hard time deciphering what it is.

[darryl]

How are you doing the matrix multiplication for the first component? The result shoud be a 1x3 vector, not 3x1. It may depend on how you're defining the vectors/matrices for exactly how the result is presented. Try doing all the matrix multiplication in one line and see if you get a single number as the result. Worse case, transpose the first result (so it's 1x3) then go as you'd planned.
Darryl

[mFidino]

So I figured out the matrix math portion, which gives us the variance in a particular site, but am getting some strange results from the rest of the process ( the upper and lower confidence intervals are both higher than the occupancy probability generated through the logistic link). Here is the process that I worked it through, using three sites as an example.

variance to standard error
SE = sqrt (x) / sqrt (y)
x = variance of a site, y = number of sites = 3

Calculate error
error = qnorm(0.975) * x / sqrt(y)
qnorm is the quantile function for a normal distribution, x = SE of a site, y = number of sites = 3

Upper and lower Confidence intervals

lowerConf = mean(x) - error
upperConf = mean(x) + error
x = covariate data

back transform

lowerProbability = 1 / (1 + e^-x)
x = lowerConf for a site

upperProbability = 1 / (1 + e^-x)
x = upperConf for a site

Results from 3 sites worth of data

mean3sites = 0.9882449
var3sites = [ 0.10577993, 0.12578995, 0.13579227 ]

SE = [ 0.1877764, 0.2047681, 0.2127536 ]

error = [ 0.2124851, 0.2317127, 0.2407489 ]

lowerConf = [ 0.7757598, 0.7565322, 0.7474960 ]
upperConf = [ 1.200730, 1.219958, 1.228994 ]

lowerProb = [1.460354, 1.469291, 1.473551 ]
upperProb = [1.300974, 1.295243, 1.292587 ]

I know I made a mistake somewhere in here, as both of these are > 1 and the probability of occupancy for the first few sites is around ~30%.

[darryl]

Why are you now talking about the mean of sites? My impression given and the approach described has been for an individual site. After doing the matrix multiplication you have the variance for the estimated probability with that certain combination of covariates.

Regardless, you clearly have a coding error. The formulas you've displayed are correct, but your transformed values from the logit to the probability scale clearly aren't.

[mFidino]

I figured out the errors that were in my code, as well as the errors in the way I had set it up. Thank you for the assistance!