Parameter counting

posts related to the RMark library, which may not be of general interest to users of 'classic' MARK

Parameter counting

Postby aswea » Fri Jan 08, 2010 12:37 pm

I'm having some trouble with knowing how to count the confounded parameter for the last recapture instance. My dataset consists of a series of recaptures of three stocks of fish each of which were released from the hatcheries in two release groups (six groups total).

With this dataset, how many confounded parameters (beta terms) would I have for the final recaptures using the following models?

1.Phi(time:STOCK)p(time:STOCK)
2. Phi(time:STOCK)p(time:RELEASE)
3. Phi(time:STOCK)p(time)

I think that Model 1 should be 3 confounded parameters-- one for each stock. However, $results$SINGULAR from RMark returned four values of which two were from recaptures mid-way through the sequence (near the boundary). That leaves only two to be the confounded terms. Why would that happen?

For Model 2, RMark seems to be telling me that there is just 1 confounded parameter (once I account for the two parameters near the boundary). Section 6.11 of the Gentle Intro (v7) describes parameter counting for additive models and says that the phi and p from the beta terms are identifiable by linear interpolation as long as we know just one. I have extended this to mean that I only have one confounded parameter for the model above because the phi and p can also be disentangled once I know just one parameter. Is this right?

It seems that it should be possible to estimate all but 1 of the confounded parameters for model 3 algebraically as well, but in the notes from my MARK course, I wrote down that there were three-- the single p value was confounded with the Phi estimate for each stock.

Thank you!

~Aswea
aswea
 
Posts: 27
Joined: Sat Oct 17, 2009 3:32 pm
Location: Gander NL

Postby jlaake » Fri Jan 08, 2010 1:09 pm

Just to be clear, RMark does NOT count parameters. It reports the parameter count assuming the model is full rank (all estimable) and the parameter count from MARK if they disagree. You can choose to use either one in model selection.

--jeff
jlaake
 
Posts: 1479
Joined: Fri May 12, 2006 12:50 pm
Location: Escondido, CA

Postby jlaake » Fri Jan 08, 2010 1:14 pm

Shouldn't have hit submit so quickly. Also, in RMark you can adjust the parameter count to be used for model selection, rather than use full rank or MARK's count. I've made the full rank count the option to avoid over-fitting which happens a lot with sparse data set and overly-complex models because many parameters will end up at boundaries and even though they are estimable MARK will not count them. I think it is better to underfit than overfit.
jlaake
 
Posts: 1479
Joined: Fri May 12, 2006 12:50 pm
Location: Escondido, CA

RE: Parameter counting

Postby dhewitt » Fri Jan 08, 2010 2:32 pm

Parameter counting seems like a hot topic these days...

aswea wrote:I'm having some trouble with knowing how to count the confounded parameter for the last recapture instance. My dataset consists of a series of recaptures of three stocks of fish each of which were released from the hatcheries in two release groups (six groups total).


I'll assume that you're dealing with a CJS (recaptures only) analysis here. If so, RMark's adjusted (full rank) count will always be wrong for fully time-dependent models -- it will be counting confounded pairs as two parameters. And if you have boundary estimates rattling around, the count from MARK will often be wrong for LOGIT link models. Bottom line: You have to know the number of estimable parameters a priori and be prepared to manually adjust some counts (as with adjust.parameter.count in RMark). Counts may match neither of the two sources. Model output gives you hints (often conclusive with good data sets) about the number of confounded parameters, but you should be checking it to be sure it matches what you know rather than vice versa.

aswea wrote:With this dataset, how many confounded parameters (beta terms) would I have for the final recaptures using the following models?

1.Phi(time:STOCK)p(time:STOCK)
2. Phi(time:STOCK)p(time:RELEASE)
3. Phi(time:STOCK)p(time)

I think that Model 1 should be 3 confounded parameters-- one for each stock. However, $results$SINGULAR from RMark returned four values of which two were from recaptures mid-way through the sequence (near the boundary). That leaves only two to be the confounded terms. Why would that happen?


I believe you are right that the model has 3 confounded parameters -- the Phi-p final pair for each stock. So, as an example assuming 4 occasions and no boundary estimates, you would count 15 parameters instead of 18. You count each pair as one estimable parameter for the model. But I think there is an issue of terminology here that returns below. It is better to think of the number of estimable parameters (which can be these confounded pairs) than the number of confounded parameters. Based on your note below about model 3 I am not entirely sure what you mean by a confounded parameter -- a pair or an individual real parameter.

I've never thought too deeply about the parameters that RMark reports as singular. IIRC, Jeff just labels them based on what MARK reports as singular. However, those are betas and not the reals, so interpreting them is a bit tricky (one confounded real parameter is not necessarily equivalent to one confounded beta, in contrast to your parenthetical note -- depends on the DM). Maybe Jeff/Gary can enlighten us on this point.

aswea wrote:For Model 2, RMark seems to be telling me that there is just 1 confounded parameter (once I account for the two parameters near the boundary). Section 6.11 of the Gentle Intro (v7) describes parameter counting for additive models and says that the phi and p from the beta terms are identifiable by linear interpolation as long as we know just one. I have extended this to mean that I only have one confounded parameter for the model above because the phi and p can also be disentangled once I know just one parameter. Is this right?


This is not an additive model, is it? If I understand the model you are referring to based on your syntax, there is surely more than 1 parameter involved in the confounding. But the terminology problem returns here... I'd say there are 6 estimable parameters in the final pairs -- the final Phi for each stock is bound up with the final p for each release (2 for each final Phi). Thus, you'd need 6 parameters for the last pairs and thus count 24 rather than 27 in our example with 4 occasions and no boundary estimates. In CJS time-dependent models, there is no way to disentangle the final Phi and p for groups. Additive models are a different story, and get nasty (IMO).

aswea wrote:It seems that it should be possible to estimate all but 1 of the confounded parameters for model 3 algebraically as well, but in the notes from my MARK course, I wrote down that there were three-- the single p value was confounded with the Phi estimate for each stock.


Given the terminology issue, I'm not entirely sure what you mean by "it should be possible to estimate all but 1 of the confounded parameters". But I think you have it right. Your count should exclude only one parameter compared to the full model parameter count. Thus, in our example (4 occasions and no boundary estimates), you would need to count 11 parameters. The final Phi for each stock is paired with the single final p.
dhewitt
 
Posts: 150
Joined: Tue Nov 06, 2007 12:35 pm
Location: Fairhope, AL 36532

Parameter counting

Postby aswea » Sat Jan 09, 2010 4:44 pm

Thank you for your reply! It is wonderful to have a place to ask these questions.

By confounded parameters, I am talking about the beta term for the final interval that is estimated because the survival and recapture rates are not individually identifiable. I think you got that, but I am not sure what made my words confusing for model 3.

I don't understand where the 27 came from in your explanation for model 2 (2. Phi(time:STOCK)p(time:RELEASE)) with the example of 4 occasions and no boundary estimates. Could you explain that please? For three stocks and two release groups, I get 15 parameters not counting the beta terms. Each release group has 1 parameter for each recapture interval (3+3) and each stock has 1 parameter for each recapture interval (3+3+3). I wish I could use a pen! Possibly what is confusing here is that I don't believe the p estimates are stock-specific. It should be the same two p estimates confounded with the final phi for each stock.

Østock1prelease1prelease2

Østock2prelease1prelease2

Østock3prelease1prelease2



or maybe it is more like:



Østock1prelease1

Østock1prelease2

Østock2prelease1

Østock2prelease2

Østock3prelease1

Østock3prelease2

In the first case, if you know any two parameters, then you can disintangle the remainder. In the second case, you only need to know one parameter to work out the others. I don't know what I am talking about, but have been trying to work out why I only lose one parameter for models like number two.

Here is the output for the betas for model two (2. Phi(time:STOCK)p(time:RELEASE)) for my data with eight recapture intervals (Group 1 joins the dataset at time5). $results$singular selects parameters 23 and 8 as confounded or unestimable. Parameter 8 is for recapture interval 5 and is clearly near the boundary. This leaves just one parameter to represent the confounding that occurs for the final interval. I am assuming that $results$singular doesn't know WHICH of the beta terms to identify as confounded and so it just highlights one.

estimate se lcl ucl
1 Phi:time5:STOCK1 1.184371 0.243334 0.707437 1.66E+00
2 Phi:time6:STOCK1 0.096179 0.289587 -0.47141 6.64E-01
3 Phi:time7:STOCK1 -1.45666 10.51368 -22.0635 1.92E+01
4 Phi:time1:STOCK2 0.934026 0.115463 0.707718 1.16E+00
5 Phi:time2:STOCK2 0.794112 0.133819 0.531827 1.06E+00
6 Phi:time3:STOCK2 1.271101 0.326584 0.630997 1.91E+00
7 Phi:time4:STOCK2 1.357937 0.423489 0.527899 2.19E+00
8 Phi:time5:STOCK2 21.53902 3163.808 -6179.52 6.22E+03
9 Phi:time6:STOCK2 -0.42277 0.276805 -0.96531 1.20E-01
10 Phi:time7:STOCK2 -1.29925 10.85429 -22.5737 2.00E+01
11 Phi:time2:STOCK3 1.058144 0.120443 0.822075 1.29E+00
12 Phi:time3:STOCK3 1.378174 0.308844 0.772839 1.98E+00
13 Phi:time4:STOCK3 1.208023 0.463944 0.298693 2.12E+00
14 Phi:time5:STOCK3 6.06774 39.46236 -71.2785 8.34E+01
15 Phi:time6:STOCK3 -0.43795 0.247057 -0.92218 4.63E-02
16 Phi:time7:STOCK3 -1.22113 11.04158 -22.8626 2.04E+01
17 p:time2:RELEASE1 3.792657 0.711862 2.397408 5.19E+00
18 p:time3:RELEASE1 4.11926 0.710048 2.727566 5.51E+00
19 p:time4:RELEASE1 -0.41881 0.173186 -0.75826 -7.94E-02
20 p:time5:RELEASE1 -2.06841 0.270189 -2.59798 ########
21 p:time6:RELEASE1 1.470544 0.257172 0.966487 1.97E+00
22 p:time7:RELEASE1 1.280575 0.470665 0.358072 2.20E+00
23 p:time8:RELEASE1 4.64468 895.0778 -1749.71 1.76E+03
24 p:time2:RELEASE2 2.776995 0.415952 1.961731 3.59E+00
25 p:time3:RELEASE2 3.060659 0.358402 2.358191 3.76E+00
26 p:time4:RELEASE2 -0.32466 0.165268 -0.64859 -7.39E-04
27 p:time5:RELEASE2 -1.94045 0.258588 -2.44728 ########
28 p:time6:RELEASE2 1.006423 0.220975 0.573312 1.44E+00
29 p:time7:RELEASE2 0.669387 0.42086 -0.1555 1.49E+00
30 p:time8:RELEASE2 0.99111 31.50109 -60.751 6.27E+01


Below is an excerpt from the output for the beta terms for model 1 (1.Phi(time:STOCK)p(time:STOCK)) which seems that it clearly should have three confounded final terms. Results$singular selects parameters 3, 10, 14 and 8 as confounded or unestimable. Two of these are for recapture interval 5 which are near boundaries and this leaves only two for the confounding that occurs for the final interval (time 7).

estimate se lcl ucl
1 Phi:time5:STOCK1 1.206017 0.289285 0.639019 1.773016
2 Phi:time6:STOCK1 -0.03067 0.334108 -0.68553 0.624178
3 Phi:time7:STOCK1 -0.34008 134.2145 -263.401 262.7204
4 Phi:time1:STOCK2 0.936467 0.115834 0.709433 1.163502
5 Phi:time2:STOCK2 0.772838 0.132307 0.513518 1.032159
6 Phi:time3:STOCK2 1.034983 0.253425 0.538271 1.531696
7 Phi:time4:STOCK2 1.475581 0.4428 0.607693 2.343469
8 Phi:time5:STOCK2 24.96187 0 24.96187 24.96187
9 Phi:time6:STOCK2 -0.25158 0.433382 -1.10101 0.597853
10 Phi:time7:STOCK2 -0.319 0 -0.319 -0.319
11 Phi:time2:STOCK3 1.086397 0.124617 0.842148 1.330646
12 Phi:time3:STOCK3 1.813403 0.554576 0.726434 2.900372
13 Phi:time4:STOCK3 0.927115 0.354843 0.231622 1.622608
14 Phi:time5:STOCK3 13.98284 1672.484 -3264.09 3292.051
15 Phi:time6:STOCK3 -0.54995 0.279622 -1.09801 -0.00189
16 Phi:time7:STOCK3 -0.19359 0 -0.19359 -0.19359


So, a part of this question is about the output from $results$singular and its interpretation.

Thanks very much for you help,

Aswea
aswea
 
Posts: 27
Joined: Sat Oct 17, 2009 3:32 pm
Location: Gander NL

Hmmm

Postby dhewitt » Fri Jan 15, 2010 8:14 pm

Sorry so late getting back to this...

I did say "If I understand the model you are referring to based on your syntax, ..." and clearly I did not understand! I assumed that "time:RELEASE" notation for the p part of the model meant that each stock-release group combo (6) had its own p by time. So you nailed it: "Possibly what is confusing here is that I don't believe the p estimates are stock-specific." It's always a little troublesome to name these interesting models in a simple way without missing key features.

So I'll admit that I am not 100% sure about this one, and frankly I am frustrated not to know. At first glance I thought there should be 16 counted parameters in agreement with your second list of the confounded pairs. That hardly makes sense since there are only 15 real parameters involved, max. However, the likelihood needs a value to stick in for individuals in each stock and each release group. I can confirm with a good data set that you only get to drop one parameter in the count with this kind of model (14 in our example), but I can't explain why. It is Friday, so maybe that has something to with it!

And, as I did before, I will defer (beg) to Jeff to explain the $results$singular issue with regard to betas and reals. It's a blur to me still, but I am sure there is a straightforward explanation. However,

I am assuming that $results$singular doesn't know WHICH of the beta terms to identify as confounded and so it just highlights one.


I am 99.9% sure this is wrong. IIRC, Jeff just labels the parameters that are flagged as singular in the MARK output (see the actual output file at the end of the real parameters list). Blame MARK, not RMark. :)
dhewitt
 
Posts: 150
Joined: Tue Nov 06, 2007 12:35 pm
Location: Fairhope, AL 36532

Postby aswea » Mon Jan 18, 2010 9:08 am

No apologies for being late! I have deadlines and questions so it is a relief to get answers and a surprise that anyone has time to provide them. Thank you!

Thanks for the confirmation of dropping one parameter with a good data set with a model like Phi(time:STOCK)p(time:RELEASE). With this, I can proceed with the most looming deadline.

When I said that $results$singular doesn’t know which of the beta terms to identify as confounded, I meant that there are several parameters for the final recapture occasion with a model such as Phi(time:STOCK)p(time:RELEASE) and yet only one of them can be identified through singular. It always chooses a survival parameter (not a recapture parameter), but not always the first one.

Yes, if anyone else can shed light on the relationship between the output of $results$singular and the betas and reals, I would be very interested. I do understand that this comes right out of MARK.

~Aswea
aswea
 
Posts: 27
Joined: Sat Oct 17, 2009 3:32 pm
Location: Gander NL

Parameter counting

Postby cschwarz@stat.sfu.ca » Wed Jan 20, 2010 6:47 pm

As dhewitt noted the use of "confounded parameters" is very confusing and you should really be thinking of estimable parameters.

Here is my take on the problem. Assume 3 stocks each of which has 2 releases.

Model 1: phi(time:stock) p(time:stock).

Let us look at the final survival * final capture terms that would appear in the likelihoods. These take the form
phi(s1,k-1) p(s1, k)
phi(s2,k-1) p(s2, k)
phi(s3,k-1) p(s3, k)
where phi(s1,k-1) is the survival term for stock 1 between sampling occasions k-1 and k; p(s1, k) is the capture rate for stock s1 at time k.

The likelihood only has 3 expressions and therefore ALL 6 parameters (3 phi and 3 p) are "confounded", i.e. cannot be estimated separately. It doesn't make sense to say that the p's are confounded but the phi's are not or vice versa. To say that there are 3 confounded parameters really doesn't have an interpretation other than perhaps saying only 3 functions of parameters can be estimated? I think the only sensible thing to say here is that there are 6 "raw" or "fundamental" parameters that cannot all be separately estimated.

So there are 6 fundamental parameters, but only 3 estimable (functions of) parameters being the 3 products above. The number of estimable parameters for these 6 end of experiment "raw" parameters is 3. When you do the parameter counting, add 3 for the final sets at the end of the experiment.


Model 2: Phi(time:stock) p(time:release). There are a total of 9 raw parameters at the end of the experiment that appear in the form:

phi(s1, k-1) p(s1-r1, k)
phi(s1, k-1) p(s1-r2, k)
phi(s2, k-1) p(s2-r1, k)
phi(s2, k-1) p(s2-r2, k)
phi(s3, k-1) p(s3-r1, k)
phi(s3, k-1) p(s3-r2, k)
where p(s1-r1,k) is the recapture rate of release 1 of stock 1 at time k, etc.

All 9 parameters are confounded and cannot be individually estimated, but 6 estimable functions of parameters can be found as follows:
(1) phi(s1, k-1) p(s1-r1, k) giving the first product above
(2) p(s1-r2)/p(s1-r1) because now [phi(s1, k-1) p(s1-r1, k)] *[p(s1-r2)/p(s1-r1)] gives the second term.
A similar patterns for stocks 2 and 3 gives estimable parameters 4...6.

Hence there are 6 estimate (functions) of the 9 raw parameters and so the parameter count for the end of study terms is only 6 -- 3 terms are lost.


Model 3: phi(time:stock) p(time). There are a total of 4 raw parameters (3 phi's and 1 p) for the end of study terms:
phi(s1, k-1) p(*,k)
phi(s2, k-1) p(*,k)
phi(s3, k-1) p(*,k)

But you can get the same 3 terms above with 3 estimable parameters being
phi(s1, k-1) p(*,k) which gives the first term above
phi(s2, k-1)/phi(s1,k-1) and the second term above can be derived as : [phi(s1, k-1) p(*,k)] * [phi(s2, k-1)/phi(s1,k-1)]
phi(s3, k-1)/phi(s1,k-1) and the third term above can be derived as: [phi(s1, k-1) p(*,k)] * [phi(s3, k-1)/phi(s1,k-1)]

So rather than counting 4 raw parameters for the end of study terms, there are only 3 estimable parameters.

In each of the cases above, the estimable functions are NOT unique and there are other functions of parameters that would give the same results, but the number of estimable functions of parameters is the same.



So how does Mark determine the number of estimable parameters. It uses a singular value decomposition of the information matrix (the inverse of the variance-covariance matrix) and looks for the number of singular values (related to eigenvalues that are close to 0). If all parameters are estimable, then the information matrix will be full rank and there will be no singular values. If the information matrix is less than full rank, this indicates that some rows (corresponding to the second derivatives of the log-likelihood) are proportional to each other and the matrix cannot be inverted (is singular). I believe this will only detect simple confounding as seen in the expressions above but won't detect "weird" confounding, but don't quote me on this. So if there is one singular value, there is one row that is a linear combination of other rows, i.e. a redundancy has been detected, but you can't associate the redundancy with a single parameter in an unambiguous way. For example, for the term phi(s1,k-1) p(s1, k-1) which of the parameters is "redundant"?

The determination of singular values is done numerically. Typically the singular values are sorted and any small singular value that is less than a very small multiple (e.g. 10**(-6)) of the largest singular value is declared to be singular. This can miss some small singular values and hence give rise to a "mis-count" of redundancies. For example, suppose the singular values were 1,000,000 1.1 .9 and the "ratio" used to determine the singular values was 10**(-6). Then the value 1.1 would NOT be declared as a "zero" singular value, but clearly is!

Now going back to the list of beta estimates that was also presented earlier, it is true that RMark only highlighted some rows, but if you examine the se for the terms, you will see that the se of the raw parameters that are involved in the confounding relationship are either 0 or enormous. Both are indications of problems regardless of what RMark flags.

These "weird" se will "identify" which parameters are involved in confounding relationships, but do NOT tell you how many estimate functions must be counted. There is no easy, automatic way way to do this (but see the papers by Catchpole and Gimenez for symbolic and numerical ways to do this respectively), and you need to write out the final terms of the likelihood and use "experience", i.e. with parameter sets as in Model 2 and 3, one estimable term is usually a ratio, to figure out what can be estimated. The first 1000 times doing this is the hardest.

So... the moral of the story is
- don't rely on automatic counting of parameters but do a systematic count based on how the PIM or DESIGN matrix is constructed.
- write out the likelihood terms near the end of the study for CJS (and related models) and similar terms near the start of the study for JS (and related models)
- use "experience" to try and figure out how many estimable functions should exist. For most CJS (and related models) these will be product terms and ratios as noted above. With experimence, you will be able to spot these fairly easily.


Carl Schwarz.
cschwarz@stat.sfu.ca
 
Posts: 43
Joined: Mon Jun 09, 2003 1:59 pm
Location: Simon Fraser University

Re: Parameter counting

Postby cooch » Wed Jan 20, 2010 7:30 pm

cschwarz@stat.sfu.ca wrote:
So how does Mark determine the number of estimable parameters. It uses a singular value decomposition of the information matrix (the inverse of the variance-covariance matrix) and looks for the number of singular values (related to eigenvalues that are close to 0). If all parameters are estimable, then the information matrix will be full rank and there will be no singular values. If the information matrix is less than full rank, this indicates that some rows (corresponding to the second derivatives of the log-likelihood) are proportional to each other and the matrix cannot be inverted (is singular). I believe this will only detect simple confounding as seen in the expressions above but won't detect "weird" confounding, but don't quote me on this. So if there is one singular value, there is one row that is a linear combination of other rows, i.e. a redundancy has been detected, but you can't associate the redundancy with a single parameter in an unambiguous way. For example, for the term phi(s1,k-1) p(s1, k-1) which of the parameters is "redundant"?

The determination of singular values is done numerically. Typically the singular values are sorted and any small singular value that is less than a very small multiple (e.g. 10**(-6)) of the largest singular value is declared to be singular. This can miss some small singular values and hence give rise to a "mis-count" of redundancies. For example, suppose the singular values were 1,000,000 1.1 .9 and the "ratio" used to determine the singular values was 10**(-6). Then the value 1.1 would NOT be declared as a "zero" singular value, but clearly is!




See - sidebar - , Chapter 4, p. 66 on. Much of what Carl summarizes (about how MARK counts parameters) is presented there in some detail. There is a *lot* of theory on estimation of parameter redundancy, and some new work recently submitted for publication which holds much promise for identifying nonidentifiable parameters regardless of underlying probability structure of the model. Most of those methods make use of high-level maths, and symbolic computation, and aren't *easily* translated to numerical routines (MARK is numerical - and does not do symbolic calculations. Some software - e.g., M-SURGE - can do more with parsing out nonidentifiable parameters, because it was built using MATLAB, which uses the MAPLE symbolic kernel for some of its routines).
cooch
 
Posts: 1652
Joined: Thu May 15, 2003 4:11 pm
Location: Cornell University

Re: Parameter counting

Postby dhewitt » Wed Jan 20, 2010 8:41 pm

Thanks a ton to Carl for weighing in on this. Models 1 and 3 are straightforward, but I still have a couple comments about Model 2...

cschwarz@stat.sfu.ca wrote:Model 2: Phi(time:stock) p(time:release). There are a total of 9 raw parameters at the end of the experiment that appear in the form:

phi(s1, k-1) p(s1-r1, k)
phi(s1, k-1) p(s1-r2, k)
phi(s2, k-1) p(s2-r1, k)
phi(s2, k-1) p(s2-r2, k)
phi(s3, k-1) p(s3-r1, k)
phi(s3, k-1) p(s3-r2, k)
where p(s1-r1,k) is the recapture rate of release 1 of stock 1 at time k, etc.


I wonder if Carl was interpreting the model name as I did originally, and incorrectly according to Aswea. I sure hope so, or I am hopelessly confused. Aswea intended the p's to be release-specific but constrained to be the same across stocks. So, there are only 2 raw p's per recapture occasion. With Carl's notation, I believe there can only be 6 raw parameters in the end because, e.g., p(s1-r1, k) is the same as p(s2-r1, k). I think this changes the game of parameter counting in this case and I haven't figured it out, even after digesting Carl's post.

cschwarz@stat.sfu.ca wrote:All 9 parameters are confounded and cannot be individually estimated, but 6 estimable functions of parameters can be found as follows:
(1) phi(s1, k-1) p(s1-r1, k) giving the first product above
(2) p(s1-r2)/p(s1-r1) because now [phi(s1, k-1) p(s1-r1, k)] *[p(s1-r2)/p(s1-r1)] gives the second term.
A similar patterns for stocks 2 and 3 gives estimable parameters 4...6.


Regardless of the previous issue, did you mean estimable parameters 3...6 here?

cschwarz@stat.sfu.ca wrote:Hence there are 6 estimate (functions) of the 9 raw parameters and so the parameter count for the end of study terms is only 6 -- 3 terms are lost.


cschwarz@stat.sfu.ca wrote:So how does Mark determine the number of estimable parameters. It uses a singular value decomposition of the information matrix (the inverse of the variance-covariance matrix) and looks for the number of singular values (related to eigenvalues that are close to 0). If all parameters are estimable, then the information matrix will be full rank and there will be no singular values. If the information matrix is less than full rank, this indicates that some rows (corresponding to the second derivatives of the log-likelihood) are proportional to each other and the matrix cannot be inverted (is singular). I believe this will only detect simple confounding as seen in the expressions above but won't detect "weird" confounding, but don't quote me on this. So if there is one singular value, there is one row that is a linear combination of other rows, i.e. a redundancy has been detected, but you can't associate the redundancy with a single parameter in an unambiguous way. For example, for the term phi(s1,k-1) p(s1, k-1) which of the parameters is "redundant"?

The determination of singular values is done numerically. Typically the singular values are sorted and any small singular value that is less than a very small multiple (e.g. 10**(-6)) of the largest singular value is declared to be singular. This can miss some small singular values and hence give rise to a "mis-count" of redundancies. For example, suppose the singular values were 1,000,000 1.1 .9 and the "ratio" used to determine the singular values was 10**(-6). Then the value 1.1 would NOT be declared as a "zero" singular value, but clearly is!

Now going back to the list of beta estimates that was also presented earlier, it is true that RMark only highlighted some rows, but if you examine the se for the terms, you will see that the se of the raw parameters that are involved in the confounding relationship are either 0 or enormous. Both are indications of problems regardless of what RMark flags.

These "weird" se will "identify" which parameters are involved in confounding relationships, but do NOT tell you how many estimate functions must be counted. There is no easy, automatic way way to do this (but see the papers by Catchpole and Gimenez for symbolic and numerical ways to do this respectively), and you need to write out the final terms of the likelihood and use "experience", i.e. with parameter sets as in Model 2 and 3, one estimable term is usually a ratio, to figure out what can be estimated. The first 1000 times doing this is the hardest.

So... the moral of the story is
- don't rely on automatic counting of parameters but do a systematic count based on how the PIM or DESIGN matrix is constructed.
- write out the likelihood terms near the end of the study for CJS (and related models) and similar terms near the start of the study for JS (and related models)
- use "experience" to try and figure out how many estimable functions should exist. For most CJS (and related models) these will be product terms and ratios as noted above. With experience, you will be able to spot these fairly easily.


Great explanation. From an experiential perspective of looking at results from models with complex design matrices (but only 673 so far!), I can support the notion that some "weird" confounding is never identified by MARK (and thus repeated in output by RMark). Also, sometimes it is not clear at all to me how the individual parameters are linked to the identified redundancies, so I would agree with that too! That was the reason for my caution to Aswea originally about the singular betas. I agree that looking at the real parameter SEs is often a giveaway, and should *reinforce* what you already suspected based on a priori logic about the design.

Thanks again to Carl and I remain curious about Model 2.
dhewitt
 
Posts: 150
Joined: Tue Nov 06, 2007 12:35 pm
Location: Fairhope, AL 36532

Next

Return to RMark

Who is online

Users browsing this forum: No registered users and 1 guest

cron