Thank you all for your help... below is the code that worked as I intended:

J <- 48
N <- 1000
S <- 0.9
PsiADD <- 0.6
PsiAA <- 0.4
PsiaA <- 0.3
p <- 0.5
c <- p

y <- matrix(0,N,J)
y[,1]="A"
dtp <- rep(c(TRUE, TRUE, FALSE, FALSE), 12)[seq_len(J)]

for (q in 1:N){
  (observable <- TRUE)
  for(j in which(!dtp)){
    if(j>1)
      if(y[q, j-1] %in% c("d", "D")) break
{
  if(j %% 2){
    if (runif(1,0,1)<=S) {
      if (observable) {
        observable <- (runif(1,0,1)<=PsiAA)
      } else {
        observable <- (runif(1,0,1)<=PsiaA)
      }
      if (observable) {
        if (runif(1,0,1) <= p) y[q,j] <- "A"
      }
    } else {
      y[q,j] <- ifelse(runif(1,0,1) <= PsiADD,"D","d")
      break
    }
  } else {
    if(observable){
      if(runif(1,0,1) <= c) y[q,j] <- "A"
    }
  }
}
  }
}

for (q in 1:N) {
  for (j in which(dtp)) {
    if(j>2)
      if(y[q, j-2] %in% c("d", "D")) break
{
  if (j %% 2) {
    if (observable) {
      observable <- runif(1, 0, 1) <= PsiAA
    } else {
      observable <- runif(1, 0, 1) <= PsiaA
    }
    if (observable) {
      if (runif(1, 0, 1) <= p) y[q, j] <- "A"
    }
  } else {
    if (observable) {
      if (runif(1, 0, 1) <= c) y[q, j] <- "A"
    }
  }
}
  }
}

On Sat, Aug 11, 2012 at 12:15 PM, Rui Barradas <ruipbarradas@sapo.pt> wrote:

> Hello,
>
> I'm not sure it works, but try the following.
>
>
> for(j in which(dtp)){
>   for (q in 1:N){
>     if(y[q, j] %in% c("d", "D")) break
>     [...etc...]
>
> and in the other loop the same,
>
>
> for (j in which(!dtp)) {
>   for (q in 1:N) {
>     if(y[q, j] %in% c("d", "D")) break
>     [...etc...]
>
>
> Em 10-08-2012 20:42, Claudia Penaloza escreveu:
>
>> This is what my code looks like now. However, there is one thing I
>> can't/don't know how to fix.
>> I can't get it to be "once dead always dead", i.e., in any given row,
>> after
>> a "D" or a "d" there should be only zeros.
>> I've tried applying a flag to break the loop if dead but I can't get it to
>> work.
>> Could you please help me with this?
>>
>> Thank you for your time,
>> Claudia
>>
>> J = 24
>> N = 10
>> S = .9
>> PsiADd = 0.4
>> PsiAA = .4
>> PsiaA = .3
>> p = .5
>> c = p
>>
>> y <- matrix(0,N,J)
>> y[,1]="A"
>> dtp <- rep(c(TRUE, TRUE, FALSE, FALSE), 6)[seq_len(J)]
>> for(j in which(dtp)){
>>    for (q in 1:N){
>>      (observable <- TRUE)
>>      if(j %% 2){
>>        if (runif(1,0,1)<=S) {
>>          if (observable) {
>>            observable <- (runif(1,0,1)<=PsiAA)
>>          } else {
>>            observable <- (runif(1,0,1)<=PsiaA)
>>          }
>>          if (observable) {
>>            if (runif(1,0,1) <= p) y[q,j] <- "A"
>>          }
>>        } else {
>>          y[q,j] <- ifelse(runif(1,0,1) <= PsiADd,"d","D")
>>          break
>>        }
>>      } else {
>>            if(observable){
>>             if(runif(1,0,1) <= c) y[q,j] <- "A"
>>            }
>>          }
>>        }
>>      }
>>
>> for (j in which(!dtp)) {
>>    for (q in 1:N) {
>>      if (j %% 2) {
>>        if (observable) {
>>          observable <- runif(1, 0, 1) <= PsiAA
>>        } else {
>>          observable <- runif(1, 0, 1) <= PsiaA
>>        }
>>        if (observable) {
>>          if (runif(1, 0, 1) <= p) y[q, j] <- "A"
>>        }
>>      } else {
>>        if (observable) {
>>          if (runif(1, 0, 1) <= c) y[q, j] <- "A"
>>        }
>>      }
>>    }
>> }
>> On Wed, Aug 8, 2012 at 2:04 PM, Claudia Penaloza
>> <claudiapenaloza@gmail.com>**wrote:
>>
>>  Answers inserted in BLUE below
>>>
>>> On Thu, Aug 2, 2012 at 5:34 PM, Claudia Penaloza <
>>> claudiapenaloza@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 <dwinsemius@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]]

