R/obj8_fitt_rgn.R

In nlr: Nonlinear Regression Modelling using Robust Methods

#+#################################################################################+|
#**       nl.fitt.rgn: Robust generalized nonlinear fitted method.                   \
#**           for storing fitted vlues                                                \
#**				Inherit from nl.fitt with all Methods.                                  /
#**       generalized form is (y-f)' V (y-f), where v(ei) = sigma^2 * V, is variance |
#**                      covariance matrix of residuals.                             |
#**			Extra Slots                                                               /
#**             vmat:    matrix part of var cov matrix of residuals.               |
#**             rmat:    choleskey decomposition of V=U'U, R = (U')^-1, this       |
#**                          choleskey decomposition is used to transfer both side  \
#**                          of model by  R*y = R*f + e, which is const varian      /
#**                                                                  --------------|
#+--------------------------------------------------+
#|      class:  nl.fitt.rob                         |
#+--------------------------------------------------+
setClass("nl.fitt.rgn",representation(
									"nl.fitt.rob",
									vm =          "matrix",
									rm =          "matrix",
									hetro =       "nl.fitt.roborNULL",
									autcorr =     "listorNULL",      ## the list can be assigned 
																	## directly to a time series fitt
									autpar =      "listorNULL",      ##  looks like repiting again some
																	## part of above time series fitt, but ok 
									gresponse =   "vectororMatrix",
									gpredictor =  "vectororMatrix"
									)
         
		)
###################################################
##   constructor
##       By default the RY and RF will be computed
###################################################
nl.fitt.rgn<-function(parameters,scale=NULL,correlation=NULL,form,
                      response =NULL,predictor,
						curvature =NULL,history =NULL,method=NULL,data,sourcefnc=NULL,Fault=Fault(),others=NULL,
						htheta=NULL,rho=NULL,ri=NULL,curvrob = NULL ,robform=NULL,	
						vm=NULL,rm=NULL,hetro=NULL,
						autcorr = NULL,autpar=NULL,
						gresponse = transformNR(response,rm),
						gpredictor = transformNR(predictor,rm)
						){
  if(is.Fault(Fault)) return(new("nl.fitt.rgn",Fault=Fault))
  result <- new("nl.fitt.rgn",
						nl.fitt.rob(
								parameters =   parameters ,
                scale =        scale,
                correlation =  correlation,
								form =         form, 
								response =     response , 
								predictor =    predictor , 
								curvature =    curvature ,
								history =      history , 
								method =       method,
								data =         data,
								sourcefnc =    sourcefnc,
								Fault =        Fault,
								others =       others,
								htheta =       htheta,
								rho =          rho,
								ri =           ri,
								curvrob =      curvrob,
								robform =      robform),
							vm =           vm,
							rm =           rm,
							hetro  =       hetro,
							autcorr =      autcorr,
							autpar =       autpar,							
							gresponse =    gresponse,
							gpredictor =   gpredictor
				)
	return(result)
}

##################################################################
##   Method residuals, nl.fitt.rgn
##      output:  residuals and gresiduals is rm . residuals
##################################################################
setMethod(f="residuals",signature="nl.fitt.rgn",
	definition=function(object,data=NULL,...){
		if(is.null(data)){
			rs <- as.numeric(object@response)
			pr <- as.numeric(object@predictor)
			rsd <- (rs-pr)
			grsd <- object@rm %*% rsd
		}
		else{
			datalist <- as.list(data)
			datalist[names(object@form@par)] <- object@parameters[names(object@form@par)]
			fmod <- eval(object@form,datalist)
			rs <- as.numeric(fmod@response)
			pr <- as.numeric(fmod@predictor)
			rsd <- (rs-pr)
			grsd <- object@rm %*% rsd
		}
	result <- structure(.Data=rsd,"gresiduals"=grsd)
	return(result)
	}
)

###################################################
##   Method parInfer, nl.fitt.rob, parameter Infrences
##      upgrade 5/12/2012 Stockholm
###################################################

setMethod(f="parInfer",signature="nl.fitt.rgn",
	definition=function(object,confidence=.95){
		if(object$method$methodID == 10) return(parInfer.WM(object))
		if(is.null(object@scale))
		{
			Fault2 <- Fault(FL=T,FN=8,FT=,nlr::db.Fault[nlr::db.Fault$FN==8,]$FT,FF="parInfer")
			#result <- Fault2
			sigma <- mscale(residuals(object))   # nl.mscale(tmp,robfunc,...)
		}
		else  sigma <- object$scale

		if(! is.null(object$hetro)) {
			varmod <- predict(object$hetro)
			varinv <- diag( 1.0 / as.numeric(varmod) )
			yhat <- predict(object)
			grdf <- attr(yhat,"gradient")
			gv <- t(grdf) %*% varinv
			gvg <- gv %*% grdf											##  (RF)' RF= F' v-1 F
			gvinv <- indifinv( gvg ,symmetric=T)
		}
		else{
			grdf <- attr(object@predictor,"gradient")
			gvg <- t(grdf) %*% grdf											##  F' F
			gvinv <- indifinv( gvg ,symmetric=T)
		}
		.expr3 <- sum(attr(object@rho,"gradient"))
		n <- nrow(grdf)
		p <- object$form$p
		spsi <- sum(.expr3^2) 
		spsid <- sum(attr(object@rho,"hessian"))
		covmatrix <- sigma^2 * ( spsi / spsid^2 ) * 
									(n*n/(n-p)) * gvinv
		v2inv <- diag(1/sqrt(diag(covmatrix)))
		corrmat <- v2inv %*% covmatrix %*% v2inv
		parstdev = sqrt(diag(covmatrix ))				

		.expr4 <- sqrt((p+1)*qf(confidence,p+1,n-p-1))
		.expr5 <- parstdev * .expr4
		cilow <- unlist(object$parameters[names(object$form$par)]) - .expr5
		ciupp <- unlist(object$parameters[names(object$form$par)]) + .expr5
	
		result <- list(covmat=covmatrix ,corrmat = corrmat,parstdev=parstdev,CI=cbind(cilow,ciupp))	
		return(result)
	}
)

#######################################################
##   Method PI, nl.fitt, prediction Interval
##     Added to this object att 05/12/2012
#######################################################
setMethod(f="predictionI",signature="nl.fitt.rgn",
	definition=function(nlfited,confidence=.95,data=NULL){
		if(is.null(data)){
      yhat <- predict(nlfited)
			data <-nlfited@data[[nlfited@form$independent]]
			grdy <-  attr(yhat,"gradient")
			yhat <- as.numeric(yhat)
			m <- nrow(grdy)
			if(! is.null(nlfited$hetro)){
			  varmod <-  predict(nlfited$hetro)
			  gt <- as.numeric(varmod)
			  grdy <- nlfited$rm %*% grdy	
			}
		}
		else{
			yhat <- predict(nlfited,newdata=data)
			grdy <- attr(yhat,"gradient")
			yhat <- as.numeric(yhat)
			data2 <- data
			##if(is.null(data[[nlfited$hetro$form$independent]])) data2[[nlfited$hetro$form$independent]] = yhat
			if(is.null(nlfited$data[[nlfited$hetro$form$independent]])) data2[[nlfited$hetro$form$independent]] = yhat
			else data2[[nlfited$hetro$form$independent]] = data[[nlfited$form$independent]]
			if(! is.null(nlfited$hetro)){
			  varmod <- predict(nlfited$hetro,newdata=data2)
			  gt <- as.numeric(varmod)
			  m <- nrow(grdy)
			  newvm <- diag(gt)
			  newrm <- diag( (1.0/sqrt(gt)))# / nlfited$parameters$sigma)
			  grdy <- newrm %*% grdy			
			}
		}

		gpg<-parInfer(nlfited)$covmat             ## sigma embded here
		gpgi <- indifinv(gpg,stp=F)
		if(is.Fault(gpgi)) gpgi <- ginv(gpg)

		t6<-rep(0,times=m)

		for(i in 1:m){
			t5<-grdy[i,]%*%gpg			#*** 1*p
			t6[i]<-t5%*%(grdy[i,])       #*** 1*1
		}

		rf <- qt((1+confidence)/2,m-nlfited@form@p)
		if(! is.null(nlfited$hetro)){ # correct in new version for autocorrelated
		  vp <- gt+t6                              ## sigma is embded inside both of variance parts
		}
		else{
		  vp <- t6                              ## sigma is embded inside both of variance parts
		}
		svp <- sqrt(vp)
		up<-yhat+rf*svp
		lw<-yhat-rf*svp
		return(list(lowerbound=lw,upperbound=up))
	}
)
#'***********************************************************
#'*   Method atypicals, nl.fitt.rgn, parameter Infrences
#'     Added to this object att 19/6/2014
#'*     corrected at 8/2/2017
#'*  According to riaz2009 et al only studres and elliptnorm
#'*     can detect outliers. 
#'*     Other measures shown here for feature development.
#'***********************************************************
setMethod(f="atypicals",signature="nl.fitt.rgn",
          definition=function(nlfited,control=nlr.control()){
            if(is.Fault(nlfited)) return(nlfited)
            result <- NULL
            v <- attr(nlfited@gpredictor,"gradient")
            hs <- attr(nlfited@gpredictor,"hessian")
            res <-as.numeric(nlfited@gresponse) -as.numeric(nlfited@gpredictor) #residuals(nlfited)
            sigma <- nlfited@scale # old one parameters[["sigma"]]
            grtgr <- t(v) %*% v      						##   F.' * F.     P*P
            ht<-eiginv(grtgr,symmetric =T)						##   (F.' F.)^-1  P*P
            vbar <- apply(v,2,mean)								##   Mean over columns of v, 
            if(class(ht)=="Fault") return(ht)
            covmat <- var(v)										##   (C) Variance covariance matrix.
            cinv <- indifinv(covmat,symmetric =T,stp=F)
            p <- nlfited$form$p
            n <- length(res)
            wnew <- potmah <- vmat <- mahd <- rep(0,n)
            
            dataset <- nlfited$data[c(nlfited$form$independent,nlfited$form$dependent)]	
            dataset <- data.frame(dataset)
            g2<-cov.mve(dataset)
            mahdata <- mahalanobis(dataset,center=g2$center,cov=g2$cov)
            mahdata <- sqrt(mahdata)
            ctfmahd1 <- median(mahdata) + 3 * mad(mahdata)	##  Mahalanobis and ctfp for all data together
            
            R <- (mahdata <= ctfmahd1)
            xr <- as.matrix(dataset[R,])						##  X, whole data, Remaining
            xtx <- t(xr) %*% xr									##  XT X
            wninv <- eiginv(xtx)									## (XT X) ^-1
            datasetx <- nlfited$data[c(nlfited$form$independent)]
            datasetx <- data.frame(datasetx)
            g3<-cov.mve(datasetx)
            mahdatax <- mahalanobis(datasetx,center=g3$center,cov=g3$cov)
            mahdatax <- sqrt(mahdatax)
            ctfmahd2 <- median(mahdatax) + 3 * mad(mahdatax) ##  Mahalanobis for only x data
            
            ##  compute vmat: square of mahalanobis distance
            for(i in 1:nrow(v)){									##   Cllasic, center and variance are classic.
              b1<-v[i,]											##  vi', will store in row, 1*p
              b2<-b1 %*% ht										##  vi' (g'g)^-1            1*p
              vmat[i] <- b2 %*% b1								##  wii = vi' (g'g)^-1 vi   1*1
              temp <- v[i,] - vbar								##  vi - mean(v)            1*p
              temp2 <- temp %*% cinv							##  (vi - vbar) * C^-1      p*1, first is row
              temp2 <- temp2 %*% (temp)						##  (vi-vbar) C^-1 (vi-vbar) 1*1 second column, 
              ##                                first row
              mahd[i] <- sqrt(temp2)							##  result after sum is a vectors
              b1 <- as.matrix(dataset[i,])
              b2<- b1 %*% wninv									##  vi' (xr'xr)^-1            1*p
              wnew[i] <- b2 %*% t(b1)							##  vi' (xr'xr)^-1 vi         1*1
            }
            ctfmahdr <- median(mahd)+3*mad(mahd)
            studres <- res / sigma / sqrt(1-vmat)				##  studentized residual, ri / (sigm sqrt(1-wii))
            elliptnorm <- studres^2 * vmat / (1-vmat) / p	##  Influence curve, Elliptic norm (Cook d)
            hadi <- vmat / (1-vmat)								##  Hadi to assess high leverage points wii/(1-wii)
            mdhad <- median(hadi)
            ctfhadi1 <- mdhad + 2 * mad(hadi)
            ctfhadi2 <- mdhad + 3 * mad(hadi)
            potmah[R] <- wnew[R] / (1 - wnew[R])				##  Potential maha lanobis Remain = wii(-d)/(1-wii(-d))
            potmah[!R] <- wnew[!R]								##  Potential for deleted data
            
            ctfpotmah <- median(potmah) + 3 * mad(potmah)	##  
            
            #####################################
            if(control$JacobianLeverage == "classic")
              JL <- jaclev(v,hs,res)           ### cllasic for robust use
            else 
              JL <- JacobianLeverage(nlfited)					###   Jacobian Leverage
            if(class(JL)!="Fault"){
              jvmat <- diag(JL)
              jl.studres <- res / sigma / sqrt(1-jvmat)					##  studentized residual, ri / (sigm sqrt(1-wii))
              jl.elliptnorm <- jl.studres^2 * jvmat / (1-jvmat) / p	##  Influence curve, Elliptic norm (Cook d)
              jl.hadi <- jvmat / (1-jvmat)								##  Hadi to assess high leverage points wii/(1-wii)
              mdhad <- median(jl.hadi)
              jl.ctfhadi1 <- mdhad + 2 * mad(jl.hadi)
              jl.ctfhadi2 <- mdhad + 3 * mad(jl.hadi)
            }
            ################################################################### output   #########################
            result1 <- structure(.Data=list(
              vmat,
              nl.robmeas(measure=studres,cutofpoint=c(2.5,3,-2.5,-3),name="Studentised Residuals"),
              nl.robmeas(measure=sqrt(elliptnorm),cutofpoint=1,name="Elliptic Norm (Cook Dist)"),
              nl.robmeas(measure=mahd,cutofpoint=c(ctfmahdr,qchisq(.95,p)),name="Regression Mahalanobis Distance"),
              nl.robmeas(measure=mahdata,cutofpoint=ctfmahd1,name="Mahalanobis MVE, data"),
              nl.robmeas(measure=mahdatax,cutofpoint=ctfmahd2,name="Mahalanobis MVE, xs"),
              nl.robmeas(measure=hadi,cutofpoint=c(ctfhadi1,ctfhadi2),name="Hadi potential"),
              nl.robmeas(measure=potmah ,cutofpoint=ctfpotmah,name="Potential mahalanobis"),
              g2,g3
            ),
            .Names=c(
              "vmat",
              "studres",
              "cook",
              "mahd.v",
              "mahd.dt",
              "mahd.xs",
              "hadi",
              "potmah",
              "mvedta","mvex"
            )
            )
            
            if(class(JL)!="Fault"){ 
              result2 <- structure(.Data=list(
                jvmat,
                nl.robmeas(measure=jl.studres,cutofpoint=c(2.5,3,-2.5,-3),name="Jacobian Leverage Studentised Residuals"),
                nl.robmeas(measure=sqrt(jl.elliptnorm),cutofpoint=1,name="Jacobian Leverage Elliptic Norm (Cook Dist)"),
                nl.robmeas(measure=jl.hadi,cutofpoint=c(jl.ctfhadi1,jl.ctfhadi2),name="Jacobian Leverage Hadi potential")
              ),
              .Names=c(
                "jl.vmat",
                "jl.studres",
                "jl.cook",
                "jl.hadi"
              )
              )
              result1[names(result2)] <- result2
            }
            
            return(result1)
          }
)


#+#################################################################################+
#|                                                                                 |
#|                   End of the object 'nl.fitt.gn'                                |
#|                                                                                 |
#|                              Sep 2009                                           |
#|                                                                                 |
#|                    Hossein Riazoshams, UPM, INSPEM                              |
#|                                                                                 |
#+#################################################################################+