www.phidot.org

Hello,

I was trying to calculate the mean of some of my Phi estimates from a CJS model using variance components, and while I read through the instructions time and time again, I cannot for the life of me figure out why my code isn't working. I have tried fixing the upper limit, as the error suggests, and it still isn't working. Any feedback would be amazing! Here is my code:

Code: Select all

global$results$AICc
zz=get.real(global,"Phi",vcv=TRUE)
z=zz$estimates$estimate[c(1,3,5,7,9,11,12)]
vcv=zz$vcv.real[c(1,3,5,7,9,11,12), c(1,3,5,7,9,11,12)]
varc=var.components(z,design=matrix(rep(1,length(z)),ncol=1),vcv, alpha = 0.05, upper = 10 * max(vcv))


Everything looks great up until the last row, where I get the error:

Error in var.components(z, design = matrix(rep(1, length(z)), ncol = 1), :
root still could not be found with increased limits; try increasing upper limit by more than a factor of 2

I have double checked that the vcv matrix "vcv" has the same number of rows as the vector containing the real estimates "z," and that my "design" matrix has the same number of rows as well. I followed the example given with the dipper data, which I can get to work in R, but can't get my own data to run.

Help! Thank you so much!
- Rebekah

Use the save function to save zz as a .rda file and send to me and I'll look at it.

Clearly I need to make the code a little smarter and a little more informative. It is possible that there is no root for the equation and that is what is occurring here. How was the model created that produced these estimates? The std error of each estimate is greater than the std deviation of the 7 estimates. The process variance is excess variability above and beyond sampling error (std error) and clearly there is none here. Also some of the estimates (1,11) and (3,12) are highly correlated >0.99. Makes me wonder how these estimates were generated. I tried with var.components.reml and it worked but I don't believe the results because the mean is greater than all of the individual estimates.


Std errors of the estimates

Code: Select all

> sqrt(diag(vcv))
         1          3          5          7          9         11         12
0.05957368 0.05104593 0.05640561 0.05985594 0.05125714 0.06068666 0.05154704


Std deviation of the estimates

Code: Select all

> sqrt(var(z))
[1] 0.03647524


#Correlation matrix

Code: Select all

> vcv/sqrt(diag(vcv))%*%t(sqrt(diag(vcv)))
           1         3         5         7         9        11        12
1  1.0000000 0.8897247 0.5216363 0.4227635 0.9072513 0.9995414 0.9232475
3  0.8897247 1.0000000 0.8535817 0.7898394 0.9992017 0.8754928 0.9968251
5  0.5216363 0.8535817 1.0000000 0.9937038 0.8320882 0.4955610 0.8093915
7  0.4227635 0.7898394 0.9937038 1.0000000 0.7647073 0.3951265 0.7384981
9  0.9072513 0.9992017 0.8320882 0.7647073 1.0000000 0.8940988 0.9992102
11 0.9995414 0.8754928 0.4955610 0.3951265 0.8940988 1.0000000 0.9111894
12 0.9232475 0.9968251 0.8093915 0.7384981 0.9992102 0.9111894 1.0000000


Hey Jeff,

Thank you for looking at my code. I generated these estimates by running a global model (in this case, the highest ranked model with the most parameters, which is not my top model) because I was informed that the global model often captures the most variance. My data is set up in a way to where I am estimating summer and winter survival for both males and females, thus I'm wanting to calculate a single mean phi estimate for summer male, summer female, winter male, and winter female, so that I can compare them in Program CONTRAST.

I ran the same var.components code for a simpler model that doesn't appear to suffer from any convergence issues (Fletcher's c-hat is close to 1, betas and SEs appear to be in order and not at boundaries), and I get the same error. I don't know what I'm doing wrong . I can send you my other code if you would like to look at it, but I wasn't intending on taking up too much of your time.

Thanks again,
Rebekah

Not sure what I can add. You can use export.MARK to export data and results from R and then import into MARK and try it's variance components. Maybe it will work. But as I mentioned, the sampling std errors are larger than the variance of the estimates you sent so there is likely no process variance that can be estimated.

Hey Jeff,

That's a great suggestion. You've given me a lot to think about, so once again, thank you for your time!

Best,
Rebekah