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.
>
> 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?
>
> # Not sure I have gotten the model assumptions down.
> # Is D" or "d" an absorbing state such as "dead or "Dead"?
>
>
>
> 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?
>
>
> 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.)
>
> 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?
>
>
>
> }
>> }
>> }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.
>>
>>
> Regards,
>>
>> Eloi
>>
>> On 12-08-01 09:47 AM, Claudia Penaloza wrote:
>>
>>> I will accept any help you can give me, especially to free myself of
>>>
>>>> SAS-brain inefficiencies...
>>>>
>>> My supervisor knows SAS not R, which may explain my code.
>>> What I'm actually trying to do is simulate mark-recapture data with a
>>> peculiar structure.
>>> It is a multistate robust design model with dummy time periods... it
>>> will basically be a matrix with a succession of letters (for different
>>> states/age classes) and zeros that are generated following a certain
>>> set of conditions (probability statements; drawing from a random
>>> uniform distribution if an event happens).
>>> Normally, survival and transition to another state occur between
>>> primary sampling periods (in my R example, every two columns.. between
>>> [,2] and [,3]) but in the "dummy time period" model survival occurs
>>> first and then transition to another state, which is why I need my
>>> "conditions" to alternate. Additionally, the robust design has
>>> secondary sampling sessions that are within the same year, i.e.,
>>> survival is assumed = 1, which is why my columns are paired. Each
>>> state can also be in an unobservable state at any given time... all of
>>> these details complicate the coding.
>>> Below I've pasted what I have thus far... I have not debugged it yet
>>> (the second loop isn't working yet and the first loop still has some
>>> glitches).
>>> Further below is properly working code for a robust design without
>>> dummy time periods (it also lacks the dead states the dummy time
>>> period model has, which happen to be part of the glitchy-ness).
>>> Is there a better (more R-ish) way of doing this?
>>> Thank you for all your help,
>>> Claudia
>>> ##############################**##############################**####
>>> # 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)
>>> 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=1)
>>> if(j %% 2){survive=runif(1,0,1)
>>> if(survive<=S){
>>> stay=runif(1,0,1)
>>> if (observable==1){
>>> if(stay<=PsiAA){
>>> observable=1
>>> }else{observable=0}
>>> }else{
>>> return=runif(1,0,1)
>>> if(return<=PsiaA){
>>> observable=1
>>> }else{observable=0}}
>>> if(observable==1){
>>> capture=runif(1,0,1)
>>> if(capture<=p){
>>> y[q,j]="A"}
>>> }}else{
>>> deado=runif(1,0,1)
>>> if(deado<=PsiAd){
>>> y[q,j]="d"
>>> }else(y[q,j]="D")
>>> if(y[q,j]%in%c("D","d")){
>>> break
>>> }
>>> }
>>> }else{
>>> if(observable==1){
>>> recap=runif(1,0,1)
>>> if(recap<=c){
>>> y[q,j]="A"
>>> }
>>> }
>>> }
>>> }
>>> }
>>>
>>> for(j in which(!dtp)){
>>> for (q in 1:N){
>>> if(j %% 2){
>>> stay=runif(1,0,1)
>>> if(observable==1){
>>> if(stay<=PsiAA){
>>> observable=1
>>> }else{observable=0}
>>> }else{
>>> return=runif(1,0,1)
>>> if(return<=PsiaA){
>>> observable=1
>>> }else{observable=0}
>>> }
>>> if(observable==1){
>>> capture=runif(1,0,1)
>>> if(capture<=p){
>>> y[q,j]="A"}
>>> }}else {
>>> if(observable==1){
>>> recap=runif(1,0,1)
>>> if(recap<=c){
>>> y[q,j]="A"
>>> }
>>> }
>>> }
>>> }
>>> }
>>> ##############################**#############
>>> ### Robust design with markovian emigration
>>> ##############################**#############
>>> J = 24
>>> N = 1000
>>> S = .9
>>> PsiOO = .4
>>> PsiUO = .3
>>> p = .5
>>> c = p
>>> y <- matrix(0,N,J)
>>> y[,1]=1
>>> for (q in 1:N){
>>> (observable=1)
>>> for(j in 2:J){
>>> if(j %% 2){surviv=runif(1,0,1)
>>> if(surviv<=S){
>>> stay=runif(1,0,1)
>>> if(observable==1){
>>> if(stay<=PsiOO){
>>> observable=1
>>> }else{observable=0}
>>> }else{
>>> return=runif(1,0,1)
>>> if(return<=PsiUO){
>>> observable=1
>>> }else{observable=0}}
>>> if(observable==1){
>>> cap=runif(1,0,1)
>>> if(cap<=p){
>>> y[q,j]=1}
>>> }else y[q,j]=0
>>> }else{break}
>>> }else{
>>> if (observable==1){
>>> recap=runif(1,0,1)
>>> if (recap<=c){
>>> y[q,j]=1}
>>> else{y[q,j]=0}
>>> }
>>> }
>>> }
>>> }
>>> On Tue, Jul 31, 2012 at 7:16 PM, David Winsemius
>>> <dwinsem...@comcast.net <mailto:dwinsem...@comcast.net**>> wrote:
>>>
>>> Claudia;
>>>
>>> The loop constructs will keep you mired in SAS-brain
>>> inefficiencies forever. Please, listen more carefully to Mercier.
>>>
>>> --
>>> David
>>>
>>>
>>>
>>> On Jul 31, 2012, at 3:54 PM, Mercier Eloi wrote:
>>>
>>> On 12-07-31 02:38 PM, Claudia Penaloza wrote:
>>>
>>> Thank you very much Rui (and Eloi for your input)... this
>>> is (the very
>>> simplified version of) what I will end up using:
>>>
>>> Could we get the extended version ? Because right now, I don't
>>> know why
>>> you need such complicated code that can be done in 1 line.
>>>
>>> J <- 10
>>> N <- 10
>>> y <- matrix(0,N,J)
>>> cols <- rep(c(TRUE, TRUE, FALSE, FALSE), 3)[seq_len(J)]
>>> for(j in which(cols)){
>>> for (q in 1:N){
>>> if(j %% 2){
>>> y[q,j]=1
>>> }else{y[q,j]=2}
>>> }
>>> }
>>> for(j in which(!cols)){
>>> for (q in 1:N){
>>> if(j %% 2){
>>> y[q,j]="A"
>>> }else{y[q,j]="B"}
>>> }
>>> }
>>>
>>> There is no need for a double loop (on N) :
>>>
>>> J <- 10
>>> N <- 10
>>> y <- matrix(0,N,J)
>>> cols <- rep(c(TRUE, TRUE, FALSE, FALSE), 3)[seq_len(J)]
>>> for(j in which(cols)){
>>> if(j %% 2){
>>> y[,j]=1
>>> }else{y[,j]=2}
>>> }
>>> for(j in which(!cols)){
>>> if(j %% 2){
>>> y[,j]="A"
>>> }else{y[,j]="B"}
>>> }
>>>
>>> if you really wants to use this code.
>>>
>>> Cheers,
>>>
>>> Eloi
>>>
>>> Cheers,
>>> Claudia
>>> On Mon, Jul 30, 2012 at 5:38 PM, Mercier Eloi
>>> <emerc...@chibi.ubc.ca <mailto:emerc...@chibi.ubc.ca>
>>> <mailto:emerc...@chibi.ubc.ca
>>>
>>> <mailto:emerc...@chibi.ubc.ca>**>> wrote:
>>>
>>> Or, assuming you only have 4 different elements :
>>>
>>> mat<- matrix(rep(c(1,2,"A", "B"),each=10),10,10, byrow=F)
>>> mat2 <- as.data.frame(mat)
>>>
>>> mat
>>>
>>> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
>>> [1,] "1" "2" "A" "B" "1" "2" "A" "B" "1" "2"
>>> [2,] "1" "2" "A" "B" "1" "2" "A" "B" "1" "2"
>>> [3,] "1" "2" "A" "B" "1" "2" "A" "B" "1" "2"
>>> [4,] "1" "2" "A" "B" "1" "2" "A" "B" "1" "2"
>>> [5,] "1" "2" "A" "B" "1" "2" "A" "B" "1" "2"
>>> [6,] "1" "2" "A" "B" "1" "2" "A" "B" "1" "2"
>>> [7,] "1" "2" "A" "B" "1" "2" "A" "B" "1" "2"
>>> [8,] "1" "2" "A" "B" "1" "2" "A" "B" "1" "2"
>>> [9,] "1" "2" "A" "B" "1" "2" "A" "B" "1" "2"
>>> [10,] "1" "2" "A" "B" "1" "2" "A" "B" "1" "2"
>>>
>>> mat2
>>> V1 V2 V3 V4 V5 V6 V7 V8 V9 V10
>>>
>>> 1 1 2 A B 1 2 A B 1 2
>>> 2 1 2 A B 1 2 A B 1 2
>>> 3 1 2 A B 1 2 A B 1 2
>>> 4 1 2 A B 1 2 A B 1 2
>>> 5 1 2 A B 1 2 A B 1 2
>>> 6 1 2 A B 1 2 A B 1 2
>>> 7 1 2 A B 1 2 A B 1 2
>>> 8 1 2 A B 1 2 A B 1 2
>>> 9 1 2 A B 1 2 A B 1 2
>>> 10 1 2 A B 1 2 A B 1 2
>>>
>>> Cheers,
>>>
>>> Eloi
>>>
>>>
>>> On 12-07-30 04:28 PM, Rui Barradas wrote:
>>>
>>> Hello,
>>>
>>> Maybe something along the lines of
>>>
>>> J <- 10
>>> cols <- rep(c(TRUE, TRUE, FALSE, FALSE), 3)[seq_len(J)]
>>> for(i in which(cols)) { do something }
>>> for(i in which(!cols)) { do something else }
>>>
>>> Hope this helps,
>>>
>>> Rui Barradas
>>>
>>> Em 31-07-2012 00 <tel:31-07-2012%2000>
>>> <tel:31-07-2012%2000>:18, Claudia Penaloza
>>>
>>> escreveu:
>>>
>>> Dear All,
>>> I would like to apply two different "for loops"
>>> to each
>>> set of four columns
>>> of a matrix (the loops here are simplifications
>>> of the
>>> actual loops I will
>>> be running which involve multiple if/else
>>> statements).
>>> I don't know how to "alternate" between the loops
>>> depending on which column
>>> is "running through the loop" at the time.
>>> ## Set up matrix
>>> J <- 10
>>> N <- 10
>>> y <- matrix(0,N,J)
>>> ## Apply this loop to the first two of every
>>> four columns
>>> ([,1:2],
>>> [,5:6],[,9:10], etc.)
>>> for (q in 1:N){
>>> for(j in 1:J){
>>> if(j %% 2){
>>> y[q,j]=1
>>> }else{y[q,j]=2}
>>> }
>>> }
>>> ## Apply this loop to the next two of every
>>> four columns
>>> ([,3:4],[,7:8],[,11:12], etc.)
>>> for (q in 1:N){
>>> for(j in 1:J){
>>> if(j %% 2){
>>> y[q,j]="A"
>>> }else{y[q,j]="B"}
>>> }
>>> }
>>> I want to get something like this (preferably
>>> without the
>>> quotes):
>>>
>>> y
>>>
>>> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
>>> [1,] "1" "2" "A" "B" "1" "2" "A" "B"
>>> "1" "2"
>>> [2,] "1" "2" "A" "B" "1" "2" "A" "B"
>>> "1" "2"
>>> [3,] "1" "2" "A" "B" "1" "2" "A" "B"
>>> "1" "2"
>>> [4,] "1" "2" "A" "B" "1" "2" "A" "B"
>>> "1" "2"
>>> [5,] "1" "2" "A" "B" "1" "2" "A" "B"
>>> "1" "2"
>>> [6,] "1" "2" "A" "B" "1" "2" "A" "B"
>>> "1" "2"
>>> [7,] "1" "2" "A" "B" "1" "2" "A" "B"
>>> "1" "2"
>>> [8,] "1" "2" "A" "B" "1" "2" "A" "B"
>>> "1" "2"
>>> [9,] "1" "2" "A" "B" "1" "2" "A" "B"
>>> "1" "2"
>>> [10,] "1" "2" "A" "B" "1" "2" "A" "B"
>>> "1" "2"
>>>
>>> Any help greatly appreciated!
>>>
>>> Claudia
>>>
>>>
>>> --
>>> Eloi Mercier
>>> Bioinformatics PhD Student, UBC
>>> Paul Pavlidis Lab
>>> 2185 East Mall
>>> University of British Columbia
>>> Vancouver BC V6T1Z4
>>>
>>>
>>>
>>>
>>> David Winsemius, MD
>>> Alameda, CA, USA
>>>
>>>
>>>
>>
>> --
>> Eloi Mercier
>> Bioinformatics PhD Student, UBC
>> Paul Pavlidis Lab
>> 2185 East Mall
>> University of British Columbia
>> Vancouver BC V6T1Z4
>>
>>
>> [[alternative HTML version deleted]]
>>
>> ______________________________**________________
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/**listinfo/r-help<https://stat.ethz.ch/mailman/listinfo/r-help>
>> PLEASE do read the posting guide http://www.R-project.org/**
>> posting-guide.html <http://www.R-project.org/posting-guide.html>
>> and provide commented, minimal, self-contained, reproducible code.
>>
>
> David Winsemius, MD
> Alameda, CA, USA
>
>
[[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.
