R/seg.Ar.fit.r

In segmented: Regression Models with Break-Points / Change-Points Estimation (with Possibly Random Effects)

seg.Ar.fit<-function(obj, XREG, Z, PSI, opz, return.all.sol=FALSE){
#-----------------
  search.min<-function(h, psi, psi.old, XREG) {
    psi.ok<- psi*h + psi.old*(1-h)
    PSI <- matrix(psi.ok, nrow=n, ncol = length(psi.ok), byrow=TRUE)
    U <- (Z - PSI) * (Z > PSI)
    #X<-cbind(XREG, U1)
    #if (pow[1] != 1) U1 <- U1^pow[1]
    obj1 <- suppressWarnings(try(eval(call.noV), silent=TRUE))
    if (class(obj1)[1] == "try-error") obj1 <- suppressWarnings(try(eval(call.noV.noinit), silent=TRUE))
    L1 <- if (class(obj1)[1] == "try-error") (L0 + 10) else (-obj1$loglik)
    L1
  }
    
    
    # call.ok    #la call con U e V
    # call.ok.noinit
    # call.noV
    # call.noV.noinit
    
  # useExp.k=TRUE
  #   #-----------------
  #   est.k<-function(x1,y1,L0){
  #       ax<-log(x1)
  #       .x<-cbind(1,ax,ax^2)
  #       b<-drop(solve(crossprod(.x),crossprod(.x,y1)))
  #       const<-b[1]-L0
  #       DD<-sqrt(b[2]^2-4*const*b[3])
  #       kk<-exp((-b[2]+ DD) /(2*b[3]))
  #       return(round(kk))
  #       
  #       #  ff<-function(xx) b[1]+b[2]*xx + b[3]*xx^2+ L0
  #       #  a<-uniroot(ff, c(log(x[4]), 3.4))
  #   }
  #   #-----------------
  #   dpmax<-function(x,y,pow=1){
  #       #deriv pmax
  #       if(pow==1) -(x>y) #ifelse(x>y, -1, 0)
  #       else -pow*((x-y)*(x>y))^(pow-1)#-pow*pmax(x-y,0)^(pow-1)
  #   }
    #-----------
    in.psi<-function(LIM, PSI, ret.id=TRUE){
        #check if psi is inside the range
        a<-PSI[1,]<LIM[1,]
        b<-PSI[1,]>LIM[2,]
        is.ok<- !a & !b #TRUE se psi e' OK
        if(ret.id) return(is.ok)
        isOK<- all(is.ok) && all(!is.na(is.ok))
        isOK}
    #------------
    far.psi<-function(Z, PSI, id.psi.group, ret.id=TRUE, fc=.93) {
        #check if psi are far from the boundaries ..s
        #   returns TRUE, if fine.
        #id.far.ok<-sapply(unique(id.psi.group), function(.x) (table(rowSums(((Z>PSI)[,id.psi.group==.x,drop=FALSE])))>=2)[-1]) #[-1] esclude lo zero, x<psi[1] 
        #id.far.ok<-sapply(unique(id.psi.group), function(.x) (tabulate(rowSums(((Z>PSI)[,id.psi.group==.x,drop=FALSE]))+1)>=2)[-1]) #[-1] esclude lo zero, x<psi[1]
        #16/01/20:
        # se un psi assume l'estremo superiore "Z>PSI" non se ne accorge, mentre Z>=PSI, si.. Il contrario e vero con estremo inf e Z>PSI
        nSeg<-length(unique(id.psi.group))
        npsij<-tapply(id.psi.group,id.psi.group,length)
        nj<-sapply(unique(id.psi.group), function(.x) { tabulate(rowSums((Z>PSI)[,id.psi.group==.x,drop=FALSE])+1) }, simplify = FALSE)    
        ff<-id.far.ok<-vector("list",length=nSeg) 
        for(i in 1:nSeg){
            if(length(nj[[i]])!=npsij[i]+1) nj[[i]]<- tabulate(rowSums((Z>=PSI)[,id.psi.group==i,drop=FALSE])+1)
            id.ok<-(nj[[i]] >= 2)
            id.far.ok[[i]] <- id.ok[-length(id.ok)] & id.ok[-1] #esattamente uguale al numero di psi del gruppo i
            ff[[i]]<-ifelse(diff(nj[[i]])>0, 1/fc, fc)
        }
        id.far.ok<-unlist(id.far.ok)
        ff<-unlist(ff)
        if(!ret.id) {return(all(id.far.ok))
        } else {
            attr(id.far.ok,"factor") <- ff
            return(id.far.ok) 
        }
        #if(ret.id) return(id.far.ok) else return(all(id.far.ok))
    } #end far.psi
    #-----------
    adj.psi<-function(psii, LIM) {pmin(pmax(LIM[1,],psii),LIM[2,])} 
    #-----------
    
    n<-nrow(Z)#length(y)
    #min.step<-opz$min.step
    rangeZ <- apply(Z, 2, range)
    alpha<-opz$alpha
    limZ <- apply(Z, 2, quantile, names=FALSE, probs=c(alpha[1], alpha[2]))
    
    digits<-opz$digits
    pow<-opz$pow
    nomiOK<-opz$nomiOK
    toll<-opz$toll
    h <-opz$h
    #conv.psi<-opz$conv.psi 
    gap<-opz$gap
    #stop.if.error<-opz$stop.if.error
    #fix.npsi<-opz$fix.npsi
    fix.npsi<-opz$stop.if.error
    dev.new<-opz$dev0
    visual<-opz$visual
    id.psi.group<-opz$id.psi.group
    it.max<-old.it.max<-opz$it.max
    psi<-PSI[1,]
    psi<-adj.psi(psi, limZ)
    PSI<-matrix(psi,nrow=n, ncol=ncol(PSI), byrow=TRUE)
    
    fc<-opz$fc
    names(psi)<-id.psi.group
    epsilon <- 10
    dev.values<-psi.values <- NULL
    #id.psi.ok<-rep(TRUE, length(psi))
    
    it <- 0
    epsilon <- 10
    k.values<-dev.values<- NULL
    psi.values <-list()
    psi.values[[length(psi.values) + 1]] <- NA
    
    #browser()
    
    nomiU<- opz$nomiU
    nomiV<- opz$nomiV
    call.ok <- opz$call.ok   #la call con U e V
    call.ok.noinit<-call.ok
    call.ok.noinit$init<-NULL
    
    call.noV <- opz$call.noV #la call con U ONLY
    call.noV.noinit<-call.noV
    call.noV.noinit$init<-NULL
    
    toll<-opz$toll
    k<-ncol(Z)
    mio.init<-NULL
    mio.init.noV<-NULL
    if(!in.psi(limZ,PSI,FALSE))  stop("starting psi out of the range.. see 'alpha' in seg.control", call.=FALSE)
    if(!far.psi(Z,PSI,id.psi.group,FALSE)) 
      stop("psi starting values too close each other or at the boundaries. Please change them (e.g. set 'quant=TRUE' 
          in seg.control()), or decrease their number.", call. = FALSE)
    n.psi1<-ncol(Z)
    #==============================================
    U <- ((Z-PSI)*(Z>PSI)) #pmax((Z - PSI), 0)^pow[1]
    colnames(U)<-nomiU
    #if(pow[1]!=1) U<-U^pow[1]
    obj0 <- suppressWarnings(try(eval(call.noV.noinit), silent=TRUE))
    if(class(obj0)[1]=="try-error")  stop("The first fit with U variables does not work..", call.=FALSE)#obj0 <- suppressWarnings(eval(call.noV.noinit)) ##a volte con i valori iniziali arima() )non converge!! Quindi provo senza init
    L0<- -obj0$loglik
    
    
    
    if(it.max==0){
      colnames(U) <- paste("U", 1:ncol(U), sep = "")
      V <- -(Z > PSI)
      colnames(V) <- paste("V", 1:ncol(V), sep = "")
      obj <- obj0 #lm.wfit(x = cbind(XREG, U), y = y, w = w, offset = offs)
      L1 <- L0
      obj$coef <- c(obj$coef, rep(0, ncol(V)))
      #names(obj$coefficients) <- names.coef
      obj$epsilon <- epsilon
      obj$it <- it
      obj <- list(obj = obj, it = it, psi = psi, psi.values = psi.values, 
                  U = U, V = V, rangeZ = rangeZ, epsilon = epsilon, nomiOK = nomiOK, 
                  SumSquares.no.gap = L1, id.psi.group = id.psi.group, 
                  id.warn = TRUE)
      return(obj)
    }
    

    
    
    n.intDev0<-nchar(strsplit(as.character(L0),"\\.")[[1]][1])
    dev.values[length(dev.values) + 1] <- opz$dev0 #del modello iniziale (senza psi)
    dev.values[length(dev.values) + 1] <- L0 #modello con psi iniziali
    psi.values[[length(psi.values) + 1]] <- psi #psi iniziali
    #==============================================
    if (visual) {
        cat(paste("iter = ", sprintf("%2.0f",0),
                  "  llik = ", sprintf(paste("%", n.intDev0+6, ".5f",sep=""), -L0), #formatC(L1,width=8, digits=5,format="f"), #era format="fg" 
                  "  k = ", sprintf("%2.0f", NA),
                  "  n.psi = ",formatC(length(unlist(psi)),digits=0,format="f"), 
                  "  ini.psi = ",paste(formatC(unlist(psi),digits=3,format="f"), collapse="  "), #sprintf('%.2f',x)
                  sep=""), "\n")
    }
    #==============================================  
    id.warn <- FALSE
    id.psi.changed<-rep(FALSE, it.max)
    pLin<- opz$pLin
    idU <- (1:ncol(U))+pLin
    idV <- (1:ncol(U))+max(idU)
    tolOp <-if(is.null(opz$tol.opt)) seq(.001, .Machine$double.eps^0.25, l=it.max) else rep(opz$tol.opt, it.max)
    #browser()
    
    
    while (abs(epsilon) > toll) {
        it<-it+1
        #if(it==2) browser()
        n.psi0 <- n.psi1
        n.psi1 <- ncol(Z)
        if(n.psi1!=n.psi0){
            U <- ((Z-PSI)*(Z>PSI)) #pmax((Z - PSI), 0)^pow[1]
            if(pow[1]!=1) U<-U^pow[1]
            obj0 <- suppressWarnings(try(eval(call.noV), silent=TRUE))
            if(class(obj0)[1]=="try-error")  obj0 <- suppressWarnings(eval(call.noV.noinit)) ##a volte con i valori iniziali arima() )non converge!! Quindi provo senza init
            L0<- -obj0$loglik
        } 
        V <- -(Z>PSI) #dpmax(Z,PSI,pow=pow[2])# ifelse((Z > PSI), -1, 0)
        X <- cbind(XREG, U, V)
        #rownames(X) <- NULL
        #colnames(X)[(ncol(XREG) + 1):ncol(X)] <- c(paste("U", 1:ncol(U), sep = ""), paste("V", 1:ncol(V), sep = ""))
        obj <- suppressWarnings(try(eval(call.ok), silent=TRUE))
        if(class(obj)[1]=="try-error")  obj <- suppressWarnings(eval(call.ok.noinit)) ##a volte con i valori iniziali arima() )non converge!! Quindi provo senza init
        beta.c <-  obj$coef[idU] #[paste("U", 1:ncol(U), sep = "")]
        gamma.c <- obj$coef[idV] #[paste("V", 1:ncol(V), sep = "")]

        if(any(is.na(c(beta.c, gamma.c)))){
            if(fix.npsi)  {
                if(return.all.sol) return(list(dev.values, psi.values)) else stop("breakpoint estimate too close or at the boundary causing NA estimates.. too many breakpoints being estimated?", call.=FALSE)
            } else {
                id.coef.ok<-!is.na(gamma.c)
                psi<-psi[id.coef.ok]
                if(length(psi)<=0) {
                    warning(paste("All breakpoints have been removed after",it,"iterations.. returning 0"), call. = FALSE)
                    return(0)
                }
                gamma.c<-gamma.c[id.coef.ok]
                beta.c<-beta.c[id.coef.ok]
                Z<-Z[, id.coef.ok, drop=FALSE]
                rangeZ <- rangeZ[,id.coef.ok, drop=FALSE]
                limZ <- limZ[,id.coef.ok, drop=FALSE]
                nomiOK<-nomiOK[id.coef.ok] #salva i nomi delle U per i psi ammissibili
                id.psi.group<-id.psi.group[id.coef.ok]
                names(psi)<-id.psi.group
            }
        }
        psi.old<-psi
        psi <- psi.old + h*gamma.c/beta.c #+h*gamma.c/beta.c
        psi<- adj.psi(psi, limZ)
        psi<-unlist(tapply(psi, opz$id.psi.group, sort), use.names =FALSE)
        
        a<-optimize(search.min, c(0,1), psi=psi, psi.old=psi.old, XREG=XREG, tol=tolOp[it])
        k.values[length(k.values) + 1] <- use.k <- a$minimum
        L1<- a$objective
        #L1.k[length(L1.k) + 1] <- L1<- a$objective
        psi <- psi*use.k + psi.old* (1-use.k)
        psi<- adj.psi(psi, limZ)
        
        if(!is.null(digits)) psi<-round(psi, digits)
        PSI <- matrix(psi, n, ncol = length(psi), byrow=TRUE)
        #--modello con il nuovo psi
        U<-(Z-PSI)*(Z>PSI) #in seg.(g)lm.fit l'ho chiamata U..
        #if(pow[1]!=1) U<-U^pow[1]
        
        if (visual) {
            flush.console()
            cat(paste("iter = ", sprintf("%2.0f",it),
                      "  llik = ", sprintf(paste("%", n.intDev0+6, ".5f",sep=""), -L1), #formatC(L1,width=8, digits=5,format="f"), #era format="fg" 
                      "  k = ", sprintf("%2.0f", k),
                      "  n.psi = ",formatC(length(unlist(psi)),digits=0,format="f"), 
                      "  est.psi = ",paste(formatC(unlist(psi),digits=3,format="f"), collapse="  "), #sprintf('%.2f',x)
                      sep=""), "\n")
        }
        
        epsilon <- (L0 - L1)/(abs(L0) + 0.1)
        L0<-L1
        #U <-U1 #le ho sempre chiamate U, non c'e' bisogno di ri-assegnare..
        
        k.values[length(k.values)+1]<-use.k
        psi.values[[length(psi.values) + 1]] <- psi
        dev.values[length(dev.values) + 1] <- L0
        
        #Mi sa che non servono i controlli.. soprattutto se non ha fatto step-halving
        #check if i psi ottenuti sono nel range o abbastanza lontani
        id.psi.far <-far.psi(Z, PSI, id.psi.group, TRUE, fc=opz$fc)
        id.psi.in <- in.psi(limZ, PSI, TRUE)
        id.psi.ok <- id.psi.in & id.psi.far 
        
        if(!all(id.psi.ok)){
            if(fix.npsi){
#                psi<-adj.psi(psi, limZ) #within range!!!
#                id.psi.far<-far.psi(Z, PSI, id.psi.group, TRUE)
                psi<-psi*ifelse(id.psi.far,1,.9) 
                PSI <- matrix(psi, n, ncol = length(psi), byrow=TRUE)
                id.psi.changed[it]<-TRUE
                    } else {
                Z<-Z[, id.psi.ok, drop=FALSE]
                PSI<-PSI[, id.psi.ok, drop=FALSE]
                rangeZ <- rangeZ[,id.psi.ok,drop=FALSE]
                limZ <- limZ[,id.psi.ok,drop=FALSE]
                nomiOK<-nomiOK[id.psi.ok] #salva i nomi delle U per i psi ammissibili
                id.psi.group<-id.psi.group[id.psi.ok]
                psi.old<- psi.old[id.psi.ok]
                psi<- psi[id.psi.ok]  
                names(psi)<-id.psi.group
                if(ncol(PSI)<=0) {
                    warning(paste("All breakpoints have been removed after",it,"iterations.. returning 0"), call. = FALSE)
                    return(0)
                }
            }
        }
        if (it >= it.max) {
            id.warn <- TRUE
            break
        }
    } #end while_it
    ##=============================================================================
    if(id.psi.changed[length(id.psi.changed)]) warning(paste("Some psi (", (1:length(psi))[!id.psi.far],
                                                             ") changed after the last iter.",sep=""), call. = FALSE)
    if(id.warn) warning(paste("max number of iterations (", it,") attained",sep=""), call. = FALSE)
    
    attr( psi.values, "dev") <- dev.values
    attr( psi.values, "k")<- k.values
    
    #ordina i breakpoints.. 
    psi<-unlist(tapply(psi, id.psi.group, sort))
    names(psi)<-id.psi.group

    names.coef<-names(coef(obj)) #names(obj$coefficients) #obj e' quello vecchio che include U1,.. V1,...
    
    #browser()
    
    #PSI.old<-PSI
    PSI <- matrix(psi, n, ncol = length(psi), byrow = TRUE)
    
    #U e V possono essere cambiati (rimozione/ordinamento psi.. ) per cui si deve ricalcolare il tutto, altrimenti sarebbe uguale a U1 e obj1
    #if(sd(PSI-PSI.old)>0 || id.psi.changed[length(id.psi.changed)]){
    V <- (Z>PSI)
    U <- (Z-PSI)*V
    V <- -V
        
    colnames(U)<-paste("U", 1:ncol(U), sep = "")
    colnames(V)<-paste("V", 1:ncol(V), sep = "")
    #X <- cbind(XREG, U, V)
    #rownames(X) <- NULL
    obj <- suppressWarnings(try(eval(call.noV), silent=TRUE))
    L1<- -obj$loglik
    obj$coef<-c(obj$coef, rep(0,ncol(V)))
    names(obj$coef)<-names.coef
    obj$epsilon <- epsilon
    obj$it <- it
    
    obj<-list(obj=obj,it=it,psi=psi, psi.values=psi.values, U=U,V=V,rangeZ=rangeZ,
              epsilon=epsilon,nomiOK=nomiOK, SumSquares.no.gap=L1, id.psi.group=id.psi.group,id.warn=id.warn) #inserire id.psi.ok?
    return(obj)
}

# if(return.all.sol) {
        #     obj.noV <- suppressWarnings(eval(call.noV)) #, envir=mfExt
        #     #mio.init.noV<-obj.noV$coef
        #     #mio.init.noV<- c(0,obj.noV$coef[-1])
        #     dev.new1 <- -obj.noV$loglik 
        #     #dev.new1 <- sum(mylm(x = cbind(XREG, U), y = y, w = w, offs = offs)$residuals^2)
        #     }
        # 
        # beta.c<-coef(obj)[nomiU]
        # gamma.c<-coef(obj)[nomiV]