### Name: sensitivty
### Title: Data Perturbations based Sensitivity Analysis
### Aliases: perturb sensitivity
### Keywords: misc optimize

### ** Examples


# Examine the sensitivity of the GLM from Venables & Ripley (2002, p.189)
# as described in the glm module.
# 
# Perturb the two independent variables using +/- 0.025
#       (relative to the size of each observations)
# uniformly distributed noise. Dependent variable is not being modified.
# 
# Summary should show that estimated coefficients are not substantively affected by noise.

if (!is.R()){
  # workaround MASS data bug in SPLUS
  require(MASS,quietly=TRUE)
}

data(anorexia,package="MASS")
panorexia = sensitivity(anorexia, glm, Postwt ~ Prewt + Treat + offset(Prewt),
     family=gaussian, 
    ptb.R=100, ptb.ran.gen=list(PTBi,PTBus,PTBus), ptb.s=c(1,.005,.005) )
summary(panorexia)

# Use classic longley dataset. The model is numerically unstable, 
# and much more sensitive to noise.  Smaller amounts of noise tremendously
# alter some of the estimated coefficients:
# 
# In this example we are not perturbing the dependent variable (employed) or 
# the year variable. So we assign then PTBi or NULL in ptb.ran.gen )

data(longley)
plongley = sensitivity(longley,lm,Employed~.) # defaults

# Alternatively, choose specific perturbation functions
#
plongley2 = sensitivity(longley,lm,Employed~., ptb.R=100, 
    ptb.ran.gen=c(list(PTBi), replicate(5,PTBus,simplify=FALSE),list(PTBi)), ptb.s=c(1,replicate(5,.001),1))

# summarizes range
sp=summary(plongley)
print(sp)
plot(sp) # plots boxplots of the distribution of the coefficients under perturbatione

# models with anova methods can also be summarized this way
anova(plongley) 

## Not run: 
##D # plots different replications
##D plot(plongley) # plots the perturbed replications individually, pausing in between
##D 
##D # plots anova results (where applicable)
##D plot(anova(plongley))
## End(Not run)

# look in summary object to extract more ...
names(attributes(sp))

# print matrix of coefficients from all runs
coef= attr(sp,"coef.betas.m")
summary(coef)



# Example where model does not accept a dataset as an argument...

# MLE does not accept a dataset as an argument, so need to
# create a wrapper function
#
# 
if (is.R()) {
     library(stats4)
     mleD<-function(data,lld,...) {
           # construct LL function with embedded data
           f=formals(lld)
           f[1]=NULL
           ll <-function()  {
              cl=as.list(match.call())
              cl[1]=NULL
              cl$data=as.name("data")
              do.call(lld,cl)
           }
           formals(ll)=f

           # call mle
           mle(ll,...)
     }

     dat=as.data.frame(cbind( 0:10 , c(26, 17, 13, 12, 20, 5, 9, 8, 5, 4, 8) ))
                                                    
     llD<-function(data, ymax=15, xhalf=6)
         -sum(stats::dpois(data[[2]], lambda=ymax/(1+data[[1]]/xhalf), log=TRUE))
     
     ## Don't show: 
     # ignore numerical warnings for example
     ow=options(warn=-1)
     
## End Don't show
     print(summary(sensitivity(dat, mleD,llD)))
     ## Don't show: 
     options(ow)
     
## End Don't show
}

# An example of using correlated noise by supplying a matrix noise function
# Note that the function can be an anonymous function supplied in the call itself

  ## Don't show: 
  if (!is.R()){
     # SPLUS ignores quietly
     ow=options(warn=-1)
  }
  
## End Don't show
  if (require(MASS,quietly=TRUE)) {
  ## Don't show: 
  if (!is.R()){
     # SPLUS ignores quietly
     ow=options(ow)
  }
  
## End Don't show

    plongleym=sensitivity(longley,lm,Employed~.,
      ptb.rangen.ismatrix=TRUE,
      ptb.ran.gen=
      function(x,size=1){
             mvrnorm(n=dim(x)[1],mu=rep(0,dim(x)[1]),
                     Sigma=matrix(.9,nrow=dim(x)[1],ncol=dim(x)[1]))*size+x}
    )
    print(summary(plongleym))
  }
  ## Don't show: 
  if (!is.R()){
     # SPLUS ignores quietly
     ow=options(ow)
  }
  
## End Don't show
 ## Don't show: 
        #unit test
        if (is.R()){
                selfTest=get("perturbTest",envir=environment(perturb))
        } else {
                selfTest=perturbTest
        }
        if (!selfTest()) {
                warning("failed self-test")
        }  
 
## End Don't show




