www.phidot.org

Hello,

When I try to run a robust design model (15 primary periods, 60secondary periods) the following script, I get only 4 p value

Code: Select all

run.RDExample=function()
{
Psi.dot=list(formula=~1)
Psi.time=list(formula=~time)
p.dot=list(formula=~1)
p.time=list(formula=~time)
gam.dot=list(formula=~1)
gam.time=list(formula=~time)
epsilon.dot=list(formula=~1)
epsilon.time=list(formula=~time)
time_intervals=c(0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0)
RD_process=process.data(aag, model="RDOccupEG", time.intervals=time_intervals)
RD_ddl=make.design.data(RD_process)

M1<-mark(RD_process, RD_ddl, model.parameters=list(Psi=Psi.dot, Gamma=gam.dot, Epsilon=epsilon.dot,  p=p.dot), invisible=TRUE)
M2<-mark(RD_process, RD_ddl, model.parameters=list(Psi=Psi.dot, Gamma=gam.dot, Epsilon=epsilon.dot,  p=p.time), invisible=TRUE)
M3<-mark(RD_process, RD_ddl, model.parameters=list(Psi=Psi.dot, Gamma=gam.dot, Epsilon=epsilon.time,  p=p.time), invisible=TRUE)
M4<-mark(RD_process, RD_ddl, model.parameters=list(Psi=Psi.dot, Gamma=gam.time, Epsilon=epsilon.time,  p=p.time), invisible=TRUE)
M5<-mark(RD_process, RD_ddl, model.parameters=list(Psi=Psi.dot, Gamma=gam.time, Epsilon=epsilon.dot,  p=p.time), invisible=TRUE)
M6<-mark(RD_process, RD_ddl, model.parameters=list(Psi=Psi.dot, Gamma=gam.time, Epsilon=epsilon.time,  p=p.dot), invisible=TRUE)
M7<-mark(RD_process, RD_ddl, model.parameters=list(Psi=Psi.dot, Gamma=gam.time, Epsilon=epsilon.dot,  p=p.dot), invisible=TRUE)
M8<-mark(RD_process, RD_ddl, model.parameters=list(Psi=Psi.dot, Gamma=gam.dot, Epsilon=epsilon.time,  p=p.dot), invisible=TRUE)
return(collect.models())
}
robustexample<-run.RDExample()


Here is the the design matrix of M4 :

Code: Select all

INPUT ---100000000000000000000000000000000

INPUT ---010000000000000000000000000000000

INPUT ---011000000000000000000000000000000

INPUT ---010100000000000000000000000000000

INPUT ---010010000000000000000000000000000

INPUT ---010001000000000000000000000000000

INPUT ---010000100000000000000000000000000

INPUT ---010000010000000000000000000000000

INPUT ---010000001000000000000000000000000

INPUT ---010000000100000000000000000000000

INPUT ---010000000010000000000000000000000

INPUT ---010000000001000000000000000000000

INPUT ---010000000000100000000000000000000

INPUT ---010000000000010000000000000000000

INPUT ---010000000000001000000000000000000

INPUT ---000000000000000100000000000000000

INPUT ---000000000000000110000000000000000

INPUT ---000000000000000101000000000000000

INPUT ---000000000000000100100000000000000

INPUT ---000000000000000100010000000000000

INPUT ---000000000000000100001000000000000

INPUT ---000000000000000100000100000000000

INPUT ---000000000000000100000010000000000

INPUT ---000000000000000100000001000000000

INPUT ---000000000000000100000000100000000

INPUT ---000000000000000100000000010000000

INPUT ---000000000000000100000000001000000

INPUT ---000000000000000100000000000100000

INPUT ---000000000000000100000000000010000

INPUT ---000000000000000000000000000001000

INPUT ---000000000000000000000000000001100

INPUT ---000000000000000000000000000001010

INPUT ---000000000000000000000000000001001


Can anyone help me?

Buco,
If you want estimates of detection to be secondary session-dependent, then you need to create a value for p.session=list(formula=~X) where X is whatever your secondary session names are (in the RDOccupancy example they are 'occ') that should provide you with secondary session specific estimates of p. Right now you have them as 'time', which has a specific meaning in RMark. I have no idea why you are getting only 4 values, the embeded design matrix for M4 is not useful as it does not tell us anything, rather we would need to see the parameter estimates (beta's and real's) at the end to get a better feel for any issues you might be having.

So, change p.time= to p.session=, use the right variable name for each session from your data, and see if that helps.

Bret

Formulas work off the design data, so when you want to understand what they are doing look at the design data. For example, str(RD_ddl$p) will show the structure of the p design data. summary(RD_ddl$p$time) will show the values of time specifically. time in an RD model is the secondary occasion value within a primary session and for each primary session they are all initiated with the same number. In this case you have 4 secondary occasions within each of the 15 primary sesssions so when you use ~time there are 4 levels to time and thus 4 parameters. Probably not a useful model in this case. If you want one p for every primary session you could use ~session to get 15 parameters and if you want a different p for each of the 60 occasions, use ~time*session.

--jeff