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 at gmail.com>wrote:
Answers inserted in BLUE below
On Thu, Aug 2, 2012 at 5:34 PM, Claudia Penaloza <
claudiapenaloza at 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 at 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"
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.
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
}
}
}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