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
<[email protected]>wrote:
> Answers inserted in BLUE below
>
> On Thu, Aug 2, 2012 at 5:34 PM, Claudia Penaloza <
> [email protected]> 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
>> <[email protected]>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]]
______________________________________________
[email protected] 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.
