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 <ruipbarra...@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
>> <claudiapenal...@gmail.com>**wrote:
>>
>> 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.
