R/robustd.0.R

In Rcapture: Loglinear Models for Capture-Recapture Experiments

"robustd.0" <- function(X, dfreq=FALSE, vt, vm="M0", vh=list("Chao"), vtheta=2, neg=TRUE)
{
############################################################################################################################################
# Validation des arguments fournis en entree et changement de leur forme si necessaire
############################################################################################################################################

        valid.one(dfreq,"logical")
        valid.vt(vt)
        Xvalid<-valid.X(X=X, dfreq=dfreq, vt=vt)
            X <- Xvalid$X
            I <- length(vt)  ## nombre de periodes primaires
        vm <- valid.vm(vm,c("none","M0","Mt","Mh","Mth"),vt,typet=FALSE)              
        vh <- valid.vh(vh,c("Chao","Poisson","Darroch","Gamma"),vm)
        vtheta <- valid.vtheta(vtheta,vh)
        valid.one(neg,"logical")
        
        
        
        
 ############################################################################################################################################
# AJUSTEMENT DU MODeLE
############################################################################################################################################

#-------------------------------------#
# elaboration et ajustement du modele #
#-------------------------------------#

        Y <- histfreq.0(X,dfreq=dfreq,vt=vt)
        fct.call <- match.call()
        histpos <- histpos.0(vt)
        Xw <- Xomega(vt,vm,vh,vtheta,fct.call,typet=FALSE,histpos)    # deuxieme composante (celle de Beta) dans le modele loglineaire du robust design
        consp<-log(apply(histpos,1,function(x,vtn){prod(choose(vtn,x))},vtn=vt))
        Xdelta<-pmin(histpos,1)
        Zw <- Zdelta(Xdelta)     # premiere composante (celle de Alpha) dans ce meme modele
        mX. <- cbind(Zw,Xw$mat)   # on fusionne ces 2 composantes explicatives
        gammanames <- paste("gamma",1:(2*I-2),sep="")
        colnames(mX.) <- c(gammanames,Xw$coeffnames)
        dimX <- dim(mX.)[2]
                
        
        # Ajustement du modele
        anaMrd <- suppressWarnings(glm(Y~offset(consp)+ mX.,family=poisson))


        # Verification du bon ajustement du modele loglineaire
        if(!anaMrd$converged) stop("'glm' did not converge")
        if(any(is.na(anaMrd$coef))) warning("some loglinear parameter estimations cannot be evaluated")
    
#-------------------------------------------------------#
# Reajustement du modele en enlevant les gamma negatifs #
#-------------------------------------------------------#

        param<-anaMrd$coef
        ppositions <- 0
        if (neg)
        {
            # Vecteur d'indicatrices pour les parametres d'interet negatifs
            indic <- as.vector(c(0,ifelse(param[2:(2*I-1)]<0,1,0)))
            for (i in 1:I)
            {
                test.Chao <- if(is.function(vh[[i]])||is.null(vh[[i]])) FALSE else if (vh[[i]]=="Chao") TRUE else FALSE
                if(test.Chao) {
                    indic <- as.vector(c(indic,rep(0,Xw$nbcoeff[i]-vt[i]+2),ifelse(param[(2*I+sum(na.rm=TRUE,Xw$nbcoeff[1:i])-vt[i]+2):(2*I+sum(na.rm=TRUE,Xw$nbcoeff[1:i])-1)]<0,1,0)))
                } else {
                    indic <- as.vector(c(indic,rep(0,Xw$nbcoeff[i])))
                }
            }
            while(sum(na.rm=TRUE,indic)>0) # Repeter la boucle jusqu'a ce qu'aucun gamma approprie ne soit negatif
            {
                # Determination de la position du premier gamma approprie negatif
                pos <- 1
                while(indic[pos]==0) pos <- pos + 1
                ppositions <- c(ppositions,pos)
                # Retrait de la bonne colonne de mX. et reajustement du modele
                mX. <- mX.[,-(pos-sum(na.rm=TRUE,ppositions<pos))]
                anaMrd <- suppressWarnings(glm(Y~offset(consp)+mX.,family=poisson))
                # Ajout de zeros dans le vecteur des parametres loglineaires
                positions <- sort(ppositions[-1])                
                param <- rep(0,dimX+1)
                param[-positions] <- anaMrd$coef 
                # Vecteur d'indicatrices pour les parametres d'interet negatifs
                indic <- as.vector(c(0,ifelse(param[2:(2*I-1)]<0,1,0)))
                for (i in 1:I)
                {
                    test.Chao <- if(is.function(vh[[i]])||is.null(vh[[i]])) FALSE else if (vh[[i]]=="Chao") TRUE else FALSE
                    if (test.Chao) {
                        indic <- as.vector(c(indic,rep(0,Xw$nbcoeff[i]-vt[i]+2),ifelse(param[(2*I+sum(na.rm=TRUE,Xw$nbcoeff[1:i])-vt[i]+2):(2*I+sum(na.rm=TRUE,Xw$nbcoeff[1:i])-1)]<0,1,0)))
                    } else {
                        indic <- as.vector(c(indic,rep(0,Xw$nbcoeff[i])))
                    }
                }
            }
        }
        positions <- sort(ppositions[-1]) 

#------------------------------------------------------------------------#
# Ajustement d'un modele pour tester la presence d'emigration temporaire #
#------------------------------------------------------------------------#

        Inono <- length(vm[vm!="none"])
        idpemig <- 0
        for (i in 2:(I-1)) if(vm[i]!="none") idpemig <- c(idpemig,i)
        idpemig <- idpemig[-1]
        
        if(Inono==3)
        {
            anaMrd2 <- NULL
            parap2 <- NULL

            mX3<-cbind(mX.,Xdelta[,idpemig])
            anaMrd3 <- suppressWarnings(glm(Y~offset(consp)+mX3,family=poisson))
            parap3<-summary(anaMrd3)$coef[length(anaMrd3$coef),1:2]
        } else if(Inono>3)
            {
                mX2<-cbind(mX.,apply(Xdelta[,idpemig],1,sum))
                anaMrd2 <- suppressWarnings(glm(Y~offset(consp)+mX2,family=poisson))
                parap2<-summary(anaMrd2)$coef[length(anaMrd2$coef),1:2]

                mX3<-cbind(mX.,Xdelta[,idpemig])
                anaMrd3 <- suppressWarnings(glm(Y~offset(consp)+mX3,family=poisson))
                nn<-dim(summary(anaMrd3)$coef)[1]
                parap3<-summary(anaMrd3)$coef[(nn-Inono+3):nn,1:2]
            } else {
                anaMrd2 <- NULL
                parap2 <- NULL
    
                anaMrd3 <- NULL
                parap3 <- NULL
            }

#---------------------------------------#
# Formation des vecteurs des parametres #
#---------------------------------------#

        # valeur de l intercept
        interc <- param[1]


        # creation des vecteurs de parametres alpha et beta
        Alpha <-rep(0,2*I-2)
        for (i in (1:(2*I-2)))  Alpha[i] <- param[i+1]

        Beta <- rep(0,length(Xw$mat[1,]))
        for (i in (1:length(Xw$mat[1,])))  Beta[i] <- param[2*I-1+i]


        # Verification de la presence de parametres gamma negatifs si l'option "neg"=FALSE
        if(!neg)
        {
            if (any(Alpha<0)) warning("one or more gamma parameters are negative,\n",
                                      "you can set them to zero with the argument 'neg'.")
        }

#--------------------------------------------------------------#
# Matrice de variances-covariances des parametres loglineaires #
#--------------------------------------------------------------#

        if(length(positions)>0) {
             varcov <- matrix(0,dimX+1,dimX+1)
             varcov[-positions,-positions] <- summary(anaMrd)$cov.unscaled
        } else varcov <- summary(anaMrd)$cov.unscaled  

############################################################################################################################################
# Estimation des parametres demographiques
############################################################################################################################################

#--------------------------------------------#
# calcul des probabilites de capture (pstar) #
#--------------------------------------------#

        pstar <- rep(0,I)
        pstarStderr <-rep(0,I)
        varcovpstar <- matrix(rep(0,I*I),ncol=I)
        dpstar<-matrix(rep(0,I*length(param)),ncol=I)
        beta <- Beta
        for (i in (1:I))
        {
                if (vm[i]=="none") # cas du no model
                {
                        pstar[i]<- exp(beta[1]+log(2^vt[i]-1))/(1+exp(beta[1]+log(2^vt[i]-1)))
                        dpstar[(2*I + if(i>1) sum(na.rm=TRUE,Xw$nbcoeff[1:(i-1)]) else 0):(2*I-1+sum(na.rm=TRUE,Xw$nbcoeff[1:i])),i] <- exp(beta[1]+log(2^vt[i]-1))/(1+exp(beta[1]+log(2^vt[i]-1)))^2
                        beta <- beta[-1]
                } else # Tous les autres modeles
                {
                        Xpf <- Xclosedp(vt[i],vm[i],vh[[i]],vtheta[i])
                        pstar[i] <- sum(na.rm=TRUE,exp(Xpf$mat%*%beta[c(1:Xpf$nbcoeff)]))/(1+ sum(na.rm=TRUE,exp(Xpf$mat%*%beta[c(1:Xpf$nbcoeff)])))
                        dpstar[(2*I + if(i>1) sum(na.rm=TRUE,Xw$nbcoeff[1:(i-1)]) else 0):(2*I-1+sum(na.rm=TRUE,Xw$nbcoeff[1:i])),i] <- t(Xpf$mat)%*%exp(Xpf$mat%*%beta[c(1:Xpf$nbcoeff)])/(1+sum(na.rm=TRUE,exp(Xpf$mat%*%beta[c(1:Xpf$nbcoeff)])))^2
                        beta <- beta[-c(1:Xpf$nbcoeff)]
                }
        }
        varcovpstar <- t(dpstar)%*%varcov%*%dpstar
        pstarStderr <- sqrt(diag(varcovpstar))    

#---------------#
# calcul des Ui #
#---------------#

        uv <- rep(1,(I-1))
        duv<-matrix(rep(0,(I-1)*length(param)),ncol=(I-1))
        for (i in (1:(I-1)))
        {
                eAlpha <- exp(Alpha[I:(I+i-1)])
                unmpstar <- (1-pstar[I:(I-i+1)])
                uv[i] <- prod(eAlpha*unmpstar)*(1-exp(-Alpha[I+i-1]))
                duv[,i] <- prod(eAlpha*unmpstar)*c(rep(0,I+i-1),exp(-Alpha[I+i-1]),rep(0,(length(param)-I-i)))
                for (j in 1:i)
                {
                    duv[,i] <- duv[,i] + (1-exp(-Alpha[I+i-1]))*prod(eAlpha[-j]*unmpstar[-j])*((1-pstar[I-j+1])*c(rep(0,I+j-1),exp(Alpha[I+j-1]),rep(0,(length(param)-I-j)))-exp(Alpha[I+j-1])*dpstar[,I-j+1])
                }
        }
        uv <- c(1,uv)
        duv <- cbind(rep(0,length(param)),duv)
        varcovuv <- t(duv)%*%varcov%*%duv
        uvStderr <- sqrt(diag(varcovuv))    

#---------------#
# calcul des Vi #
#---------------#

        vv <- rep(1,(I-1))
        dvv<-matrix(rep(0,(I-1)*length(param)),ncol=(I-1))
        for (i in (1:(I-1)))
        {
                eAlpha <- exp(Alpha[1:i])
                unmpstar <- (1-pstar[1:i])
                vv[i] <- prod(eAlpha*unmpstar)*(1-exp(-Alpha[i]))
                dvv[,i] <- prod(eAlpha*unmpstar)*c(rep(0,i),exp(-Alpha[i]),rep(0,(length(param)-i-1)))
                for (j in 1:i)
                {
                    dvv[,i] <- dvv[,i] + (1-exp(-Alpha[i]))*prod(eAlpha[-j]*unmpstar[-j])*((1-pstar[j])*c(rep(0,j),exp(Alpha[j]),rep(0,(length(param)-j-1)))-exp(Alpha[j])*dpstar[,j])
                }
        }
        vv <- c(1,vv)
        dvv <- cbind(rep(0,length(param)),dvv)
        varcovvv <- t(dvv)%*%varcov%*%dvv
        vvStderr <- sqrt(diag(varcovvv))    

#--------------------------------------------------------------#
# calcul des probabilites de survie entre chaque periode (phi) #
#--------------------------------------------------------------#

        phi <- rep(0,(I-1))
        phistderr <-rep(0,(I-1))
        varcovphi <- matrix(rep(0,(I-1)*(I-1)),ncol=(I-1))
        dphi<-matrix(rep(0,(I-1)*length(param)),ncol=(I-1))
        
        phi[(I-1):1] <-1/(1+uv[2:I]/cumsum(uv[1:(I-1)]))   
        
        for ( i in 1:(I-1))
        {
            if(i==I-1)
            {
                dphi[,i] <- -(phi[i]^2)*duv[,I-i+1]
            } else {
                dphi[,i] <- -(phi[i]^2)*(duv[,I-i+1]*sum(na.rm=TRUE,uv[1:(I-i)])-uv[I-i+1]*rowSums(duv[,1:(I-i)]))/sum(na.rm=TRUE,uv[1:(I-i)])^2
            }
        }
        varcovphi <- t(dphi)%*%varcov%*%dphi
        phistderr <- sqrt(diag(varcovphi))    
  
#---------------------------------------#
# calcul des taille de populations (Ni) #
#---------------------------------------#

        Npop <- rep(0,I)
        NpopStderr <-rep(0,I)
        varcovtpop <- matrix(rep(0,I*I),ncol=I)
        dNpop<-matrix(rep(0,I*length(param)),ncol=I)

        Npop[1] <- exp(interc)/(prod(1-pstar)*prod(phi))
        dprodpstar <-rep(0,length(param))
        for (i in 1:I)
        {
            dprodpstar<-dprodpstar-prod(1-pstar[-i])*dpstar[,i]
        }
        dprodphi <-rep(0,length(param))
        for (i in 1:(I-1))
        {
            dprodphi<-dprodphi+prod(phi[-i])*dphi[,i]
        }
        dNpop[,1] <- (prod(1-pstar)*prod(phi)*c(exp(interc),rep(0,(length(param)-1))) - exp(interc)*(prod(phi)*dprodpstar+prod(1-pstar)*dprodphi))/ (prod(1-pstar)*prod(phi))^2
        
        Npop[2:I]<-Npop[1]*cumprod((vv[2:I]/cumsum(vv[1:(I-1)])+1)*phi)
        for (i in 2:I)
        {
            if(i==2) {
                dNpop[,i] <- phi[i-1]*(vv[i]+1)*dNpop[,i-1] + Npop[i-1]*(vv[i]+1)*dphi[,i-1] + Npop[i-1]*phi[i-1]*dvv[,i]          
            } else {
                dNpop[,i] <- phi[i-1]*(vv[i]/sum(na.rm=TRUE,vv[1:(i-1)])+1)*dNpop[,i-1] + Npop[i-1]*(vv[i]/sum(na.rm=TRUE,vv[1:(i-1)])+1)*dphi[,i-1] + Npop[i-1]*phi[i-1]*(sum(na.rm=TRUE,vv[1:(i-1)])*dvv[,i]-vv[i]*rowSums(dvv[,1:(i-1)]))/sum(na.rm=TRUE,vv[1:(i-1)])^2
            }
        }
        varcovtpop <- t(dNpop)%*%varcov%*%dNpop
        NpopStderr <- sqrt(diag(varcovtpop))    
        NpopStderr <- sqrt(pmax(NpopStderr^2-Npop,0))

#----------------------------#
# calcul des naissances (Bi) #
#----------------------------#

        B<-Npop[2:I]-Npop[1:(I-1)]*phi
        dB <- dNpop[,2:I] - t(phi*t(dNpop[,1:(I-1)])) - t(Npop[1:(I-1)]*t(dphi))
        varcovB <- t(dB)%*%varcov%*%dB
        BStderr <- sqrt(diag(varcovB))
                
#--------------------------------------------------------------------#
# Calcul du nombre total d'individus qui ont passe sur le territoire #
#--------------------------------------------------------------------#

        Ntot <- Npop[1]+sum(na.rm=TRUE,B)
        dNtot <- dNpop[,1] + rowSums(dB)    
        NtotStderr <- sqrt(max(t(dNtot)%*%varcov%*%dNtot-Ntot,0))

############################################################################################################################################
# Presentation des resultats
############################################################################################################################################

        modelfit <- matrix(c(anaMrd$deviance,anaMrd$df.residual,anaMrd$aic),nrow=1)
        dimnames(modelfit) <- list("fitted model",c("deviance","    df","      AIC"))
        
        emigfit <- cbind(c(anaMrd2$deviance,anaMrd3$deviance),c(anaMrd2$df.residual,anaMrd3$df.residual),c(anaMrd2$aic,anaMrd3$aic))
        if (Inono==3) { dimnames(emigfit) <- list(c("model with temporary emigration"),c("deviance","    df","      AIC"))
        } else if (Inono>3) { dimnames(emigfit) <- list(c("model with homogeneous temporary emigration","model with temporary emigration"),c("deviance","    df","      AIC"))}    
               
        titre.periode.emig<-rep(0,length(idpemig))
        for (i in 1:length(idpemig)){ titre.periode.emig[i]<-paste("period",idpemig[i])}
        titre.periode<-rep(0,I)
        for (i in 1:I){titre.periode[i]<-paste("period",i)}
        titre.inter.periode<-rep(0,I-1)
        for (i in 1:(I-1)){titre.inter.periode[i]<-paste("period",i,"->",i+1)}
        titre.i<-rep(0,I)
        for (i in 1:I){titre.i[i]<-paste("i =",i-1)}
             
        parap <- rbind(parap3,parap2)     
        pstar <- cbind(pstar,pstarStderr)
        phi <- cbind(phi,phistderr)
        Npop <- cbind(Npop,NpopStderr)
        B <- cbind(round(B,digits=6),round(BStderr,digits=6))
        Ntot <- cbind(Ntot,NtotStderr)
        loglinearpara <- cbind(param,sqrt(diag(varcov)))
        uv <- cbind(uv,uvStderr)
        vv <- cbind(vv,vvStderr)

        if (Inono==3) { dimnames(parap) <- list(titre.periode.emig,c("estimate","stderr")) 
        } else if (Inono>3) { dimnames(parap) <- list(c(titre.periode.emig,"homogenous"),c("estimate","stderr")) }
        dimnames(pstar)<-list(titre.periode,c("estimate","stderr"))
        dimnames(phi)<-list(titre.inter.periode,c("estimate","stderr"))
        dimnames(Npop)<-list(titre.periode,c("estimate","stderr"))
        dimnames(B)<-list(titre.inter.periode,c("estimate","stderr"))
        dimnames(Ntot)<-list("all periods",c("estimate","stderr"))
        dimnames(loglinearpara) <- list(c("intercept",gammanames,Xw$coeffnames),c("estimate","stderr"))
        dimnames(uv) <- list(titre.i,c("estimate","stderr"))
        dimnames(vv) <- list(titre.i,c("estimate","stderr"))

        models<-matrix(Xw$models,nrow=1)
        dimnames(models) <- list("model",titre.periode)

        # Matrice de variances-covariances des parameters pstar, phi, Npop, B et Ntot
        dP <- cbind(dpstar,dphi,dNpop,dB,dNtot)
        covP <- t(dP)%*%varcov%*%dP
        titre.P <- c(rep(0,4*I-2),"Ntot")
        for (i in 1:I){
            titre.P[i]<-paste("p*",i)
            titre.P[2*I-1+i] <- paste("Npop",i)
        }
         for (i in 1:(I-1)){
            titre.P[I+i]<-paste("phi",i)
            titre.P[3*I-1+i] <- paste("B",i)
        }
        dimnames(covP)<-list(titre.P,titre.P)

     
        ans<-list(n=sum(na.rm=TRUE,Y),models=models,model.fit=modelfit,emig.fit=emigfit,emig.param=parap,
                  capture.prob=pstar,survivals=phi,N=Npop,birth=B,Ntot=Ntot,
                  loglin.param=loglinearpara,u.vector=uv,v.vector=vv,cov=covP,neg=positions)
        class(ans) <- "robustd"
        ans        
}