Answers inserted in BLUE below
On Thu, Aug 2, 2012 at 5:34 PM, Claudia Penaloza
<claudiapenal...@gmail.com>wrote:
> Thank you very much for all your suggestions. I am very sorry my code is
> so crude (it gives me a headache too!), I'm very new to R programing. I
> will have to look at your suggestions/questions very carefully (I'm barely
> holding on at this point) and get back to you (Dr. Winsemius) with some
> answers.
>
> Thank you!
> Claudia
>
> On Wed, Aug 1, 2012 at 2:11 PM, David Winsemius <dwinsem...@comcast.net>wrote:
>
>>
>> On Aug 1, 2012, at 12:02 PM, Mercier Eloi wrote:
>>
>> Your code almost gave me a headache. :-/
>>>
>>
>> I had a similar reaction. However, my approach might have been to request
>> a more complete verbal description of the data structures being operated on
>> and the methods and assumptions being used. Generally it is going to be
>> easier to go from a procedural description to good R code than it is to go
>> from a SAS Data Proc ... even if it were tested and debugged... and yours
>> was not even debugged.
>>
>>
>> ##############################**##############################**####
>>> # Robust design with Markovian emigration and dummy time periods
>>> ##############################**##############################**####
>>> J = 24
>>> N = 10
>>> S = .9
>>> PsiAD = 0.06
>>> PsiAd = 0.04
>>> PsiAA = .4
>>> PsiaA = .3
>>> p = .5
>>> c = p
>>> y <- matrix(0,N,J)
>>>
>>
>> # So is 'y' a matrix of states where the N rows are levels and the J
>> columns are discrete times?
>> # Actually I decided that context suggested you were using letters to
>> represent states.
>>
> Yes, 'y' is a matrix of state, N rows are levels (independent of each
other) and J columns are discrete times.
Yes, I am using letters to represent states.
>
>> y[,1]="A"
>>>
>>
>> # So we start with N <something>'s in state "A"?
>> # It seems as though it might be the case that every row is independent
>> of the others.
>> # .. and you are performing ( replicate()-ing in R terminology) this test
>> N times
>> # states:
>> #("A" "D")
>> #("a" "d")
>> #transitions:
>> matrix( c( runif(1, 0, 1) <= 0.4, # A -> A
>> runif(1,0,1) <= 0.3, # a -> A
>> runif(1,0,1) <= 0.06. # A -> D
>> runif(1,0,1) <= 0.04 # A -> d) )
>>
>> # What is state "a"?
>> # How do you get from A to a?
>> # can you get from D or d to A or a?
>>
> Yes, we start with N "individuals" in state "A", each individual is
independent of each other.
State "a" is an unobserved (!observable) version of state "A"
(biologically, an individual that has temporarily left the sampling area
and is not available for capture)
An individual in state "A" transitions to state "a" if it is observable
and doesn't stay (stay >= AA)
"D" and "d" are dead (absorbing) states. Once in either "D" or "d", the
individual can no longer transition and is no longer "captured" (the row
should only have zeros after a "D" or a "d")
>
>
>
>> # Not sure I have gotten the model assumptions down.
>> # Is D" or "d" an absorbing state such as "dead or "Dead"?
>>
> Yes.
> dtp <- rep(c(TRUE, TRUE, FALSE, FALSE), 6)[seq_len(J)]
>>
>>
>> #Presumably the "dummy time periods"
>>
> Yes.
>
>>
>> for(j in which(dtp)){ # j will be c(1,2, 5,6, 9,10 , ..., 21,22)
>>
>>
>>
>> for (q in 1:N){ # This can almost surely become vectorized
>>
>> (observable=1)
>> if(j %% 2){survive=runif(1,0,1)
>> if(survive<=S){
>> stay=runif(1,0,1)
>> if (observable==1){
>>
>>
>> # It would help if the concept of "observable" were explained
>> # Is this something to do with the capture-recapture observables?
>>
>
Yes. All individuals start out "observable" (we need to capture them a
first time). Individuals can then stay in their current state or transition
to another one, if they are observable they can continue to be observable,
or they can become unobservable and viceversa. Transition depends on what
state you are in at any given time (A->A != A->a != a->A != a->a).
>
>>
>> if{
>> observable=1
>> }else{observable=0}
>> }else{
>>
>>
>> # After allowing for Mercier's astute observation that observable will
>> always be 1, I'm wondering if it can be replaced with just this code (and
>> not need the loop surrounding it.)
>>
>>
>>
> The individual will always "start" as observable... but as the loop
progresses through the columns in a row, it can go between being observable
and unobservable (at least that is what I wanted to code).
>
>>
>> observable <- (stay<=PsiAA)
>>
>> #--------------
>>
>> return=runif(1,0,1)
>>
>>
>> # better NOT to use the name "return" since it is a crucial R function
>> name.
>>
>>
>>
> Will change. Thank you.
>
>>
>>
>>
>> if(return<=PsiaA){
>> observable=1
>> }else{observable=0}}
>>
>>
>> #------- perhaps:
>>
>> randret <- return=runif(N,0,1)
>> observable <- randret <= PsiaA
>> # -------------------
>>
>>
>> if(observable==1){
>> capture=runif(1,0,1)
>> if(capture<=p){
>>
>>
>> #---------- perhaps:
>> randcap <- runif(N, 0,1)
>> y[ randcap <= p, j] <- "A"
>> #That would replace with "A" (within the j-column) ...
>> # only the rows in 'y' for which randcap were less than the randcap
>> threshold.
>>
>> #------------
>>
>>
>> y[q,j]="A"}
>> }}else{
>>
>>
>>
>>
>> deado=runif(1,0,1)
>> if(deado<=PsiAd){
>> y[q,j]="d"
>> }else(y[q,j]="D")
>>
>>
>> # -------Perhaps
>> deado <- runif(N, 0,1)
>> y[ , j] <- ifelse( deado<=PsiAd, "d", "D")
>> #------------------------
>>
>>
>>
>> if(y[q,j]%in%c("D","d")){
>> break
>>
>>
>> # ----------Really? I thought that condition was assured at this
>> point?
>>
>>
>>
> I thought so too but non-zero elements appeared toward the right in rows
that had already "died".
>
>> }
>> }
>> }else{
>> if(observable==1){
>> recap=runif(1,0,1)
>> if(recap<=c){
>> y[q,j]="A"
>> }
>> }
>> }
>> }
>> }
>>
>>
>>
>>
>>
>>
>> There are a lot of unnecessary tests and conditions. I tried to break
>> down the code and write the tests that have been done when assigning a
>> variable. I simplified your the first part but cannot guarantee that all
>> criteria are met.
>>
>>
>>
> I will run the three modified versions and see how things go. I hope my
answers are helpful.
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
