Nothing
#+#################################################################################+|
#** 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 |
#| |
#+#################################################################################+
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.