predictLogProb {STAR} | R Documentation |
This function is designed to select models by cross-validation.
If two models A and B managed to pass the Ogata's tests, one should split the
data set in two parts, part1 and part2 fit each model on each part to get four
fitted model objects: A1, A2, B1, B2. The chosen model should then be the one
giving the largest of: predictLogProb(A1,part2) + predictLogProb(A2,part1)
and
predictLogProb(B1,part2) + predictLogProb(B2,part1)
.
predictLogProb(object, newdata)
object |
an object inheriting from ssanova
and ssanova0 (gssanova and
gssanova0 objects are therefore suitable). |
newdata |
a data frame containing the required variables. This
data frame must be different from the one used to obtain object . |
If eta[i] is the prediction of the fitted model for element i of
newdata
the log probability is given by :
event[i] * eta[i] - log(1 + exp(eta[i]))
Where event[i]
is 0 or 1 depending on the absence or presence of
a spike at the considered time.
A binomial regression is assumed here.
A numeric, the sum over the index i
above, the log probability
of the data contained in newdata
assuming that the model contained in object
is correct.
Christophe Pouzat christophe.pouzat@gmail.com
mkGLMdf
,
quickPredict
gssanova
## Not run: data(e060824spont) summary(e060824spont[["neuron 1"]]) reportHTML(e060824spont[["neuron 1"]],filename="e060824spont_1",otherST=e060824spont[c(2)],maxiter=100) acf(diff(e060824spont[["neuron 1"]]),type="partial") DFA <- subset(mkGLMdf(e060824spont,0.004,0,59),neuron==1) DFA <- within(DFA,i1 <- isi(DFA,lag=1)) DFA <- DFA[complete.cases(DFA),] m2u1 <- mkM2U(DFA,"lN.1",0,29,seed=20061001) m2ui <- mkM2U(DFA,"i1",0,29,maxiter=200,seed=20061001) DFA <- within(DFA,e1t <- m2u1(lN.1)) DFA <- within(DFA,i1t <- m2ui(i1)) with(DFA,plot(ecdf(e1t),pch=".")) with(DFA,plot(ecdf(i1t),pch=".")) DFAts <- as.ts(apply(DFA[,c("e1t","i1t")],2,qnorm)) plot(filter(DFAts,rep(1/125,125))) system.time(GF1 <- gssanova(event ~ e1t+i1t, data=subset(DFA,time<=29),family="binomial",seed=20061001)) tt.1 <- GF1 %tt% subset(DFA,time>29) tt.1.summary <- summary(tt.1) tt.1.summary plot(tt.1.summary,which=c(1,2,4,6)) renewalTestPlot(tt.1$ppspFct()) plot(GF1,nc=1,nr=2) system.time(GF2 <- gssanova(event ~ e1t+i1t, data=subset(DFA,time>29),family="binomial",seed=20061001)) tt.2 <- GF2 %tt% subset(DFA,time<=29) tt.2.summary <- summary(tt.2) tt.2.summary plot(tt.2.summary,which=c(1,2,4,6)) renewalTestPlot(tt.2$ppspFct()) plot(GF2,nc=1,nr=2) system.time(GF3 <- gssanova(event ~ e1t*i1t, data=subset(DFA,time<=29),family="binomial",seed=20061001)) tt.3 <- GF3 %tt% subset(DFA,time>29) (tt.3.summary <- summary(tt.3)) plot(tt.3.summary,which=c(1,2,4,6)) renewalTestPlot(tt.3$ppspFct()) plot(GF3,nc=1,nr=3) system.time(GF4 <- gssanova(event ~ e1t*i1t, data=subset(DFA,time>29),family="binomial",seed=20061001)) tt.4 <- GF4 %tt% subset(DFA,time<=29) (tt.4.summary <- summary(tt.4)) plot(tt.4.summary,which=c(1,2,4,6)) renewalTestPlot(tt.4$ppspFct()) plot(GF4,nc=1,nr=3) ## Get the log probability of the data with the additive model predictLogProb(GF1,newdata=subset(DFA,time>29))+predictLogProb(GF2,newdata=subset(DFA,time<=29)) ## Get the log probability of the data with the non-additive model predictLogProb(GF3,newdata=subset(DFA,time>29))+predictLogProb(GF4,newdata=subset(DFA,time<=29)) ## The non additive model is the "best" so refit it to the whole data set system.time(GF5 <- gssanova(event ~ e1t*i1t, data=DFA,family="binomial",seed=20061001)) plot(GF5,nr=3,nc=1) ## End(Not run)