[R] Error with stepwise models selection using NADA library

Abby Spurdle @purd|e@@ @end|ng |rom gm@||@com
Tue Mar 10 10:13:42 CET 2020


I'm assuming that you're referring to MASS::stepAIC.

There's a few problems here.
Firstly, I don't think stepAIC is compatible with S4 objects.

My first strategy was to create another class.
However, that requires many new methods.
(And the stepAIC documentation doesn't say what methods are needed).

My second strategy was to wrap the cenreg function and return a modified
copy of the survreg object.
This appears to work, but there are complications using the update and
terms functions.

Also, I got the impression that the return value from calling the term
function on the NADA-based objects, is also not compatible with the stepAIC
function, possible causing the S4 error messages.

I got the following working, hopefully...

----------
library (NADA)
library (MASS)
data (TCEReg)

cenreg.auto <- function (formula, ...)
{   mfit <- cenreg (formula, ...)@survreg
    class (mfit) <- c ("cenreg.auto", class (mfit) )
    mfit$call <- match.call ()
    mfit
}

update.cenreg.auto <- function (object, formula, ..., evaluate=TRUE)
{   call <- object$call
    f0 <- eval (call [[2]])
    f0 <- update.formula (f0, formula)
    call [[2]] <- as.call (f0)
    if (evaluate)
    {    object <- eval (call)
        object$call <- call
        object
    }
    else
        call
}

terms.cenreg.auto <- function (x, ...)
{  terms <- main.formula
   attributes (terms) <- attributes (x$terms)
   terms
}

#main dataset
#(current code needs this to be *global*)
attach (TCEReg)

#simulated data
#model sometimes picks these up
#(correct or bug?)
xrand <- 50 * sample (LandUse)
yrand <- 50 * sample (Depth)

#main formula
#(current code needs this to be *global* too)
main.formula <- Cen (TCEConc, TCECen) ~ LandUse + xrand + Depth + yrand

#fit
#(current code needs *formula duplicated*)
myfit <- cenreg.auto (Cen (TCEConc, TCECen) ~ LandUse + xrand + Depth +
yrand)

bestfit.if.lucky <- stepAIC (myfit, Cen (TCEConc, TCECen) ~ .)

detach (TCEReg)

print (bestfit.if.lucky)
----------

There are some problems with this code:
(1) Both the dataset and formula need to be in the global environment.
(2) The formula needs to be typed twice.

I don't consider this good R code.
However, I'm not an expert on formula objects, and I spent more time on
this problem than planned.
Someone who is an expert on formula objects could improve this,
substantially.

The other thing to note, is that I included some terms in the input, that
contain random numbers.
They appear in the final model, more frequently than I initially expected.
However, I'm not familiar with the type regression that you're using, so
can't really say how often this should happen.


On Sat, Mar 7, 2020 at 1:21 AM Marco Gabrielli <marco.gabrielli94 using gmail.com>
wrote:

> Hi,
> I'm working with environmental data which are left-censored and I found the
> R package NADA which seems to do the job.
> After fitting a complete model, using the cenreg function,I'd like to do a
> stepwise procedure to minimize the AIC using the stepAIC function.
> Unfortunately, after fitting the model, the stepAIC does not work with
> either the cenreg object or the S3 survreg object within.
>
> Here it is an example of the commands used:
>
> library(NADA)
> MLEtot=with(df, cenreg(Cen(Y,Cen)~x1+x2+x3+x4,dist="lognormal"))
> stepAIC(MLEtot)
> "Error: $ operator not defined for this S4 class"
> stepAIC(MLEtot using survreg)
> "Start:  AIC=6.16
> asSurv(Cen(Y, Cen)) ~ x1 + x2 + x3 + x4
> Error in (function (classes, fdef, mtable)  :
>    unable to find an inherited method for function ‘asSurv’ for signature
> ‘"Surv"’"
>
> How can I do manage to make it work?
> Thanks
>
> ______________________________________________
> R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

	[[alternative HTML version deleted]]



More information about the R-help mailing list