Nothing
covariateModelSelection <- function(pen.coef=NULL, weight=1, n.full=10,
nb.model=1, covToTransform=NULL, covFix=NULL,
direction="both", paramToUse="all", steps=1000, p.max=1,
sp0=NULL, iter=1, correlation.model=NULL, eta=NULL) {
# correlation.model=NULL
# project.folder <- mlx.getProjectSettings()$directory
# sp.file <- file.path(project.folder,"IndividualParameters","simulatedIndividualParameters.txt")
# sp.df <- read.res(sp.file)
sp.df <- mlx.getSimulatedIndividualParameters()
if (is.null(sp.df$rep))
sp.df$rep <- 1
# sp.df <- sp.df %>% arrange(rep)
if (!is.null(sp0)) {
nrep0 <- max(sp0$rep)
sp.df$rep <- sp.df$rep + nrep0
dn <- setdiff(names(sp.df), names(sp0))
sp0[dn] <- sp.df[dn]
sp.df <- rbind(sp0, sp.df)
}
nrep <- max(sp.df$rep)
ind.dist <- mlx.getIndividualParameterModel()$distribution
param.names <- names(ind.dist)
# if (identical(paramToUse,"all"))
# paramToUse <- param.names
n.param <- length(param.names)
#sim.parameters <- sp.df[c("rep","id",param.names)]
# sim.parameters <- mlx.getSimulatedIndividualParameters()
cov.info <- mlx.getCovariateInformation()
cov.names <- cov.info$name
cov.types <- cov.info$type
# j.strat <- grep("stratification",cov.types)
# if (length(j.strat)>0) {
# cov.names <- cov.names[-j.strat]
# cov.types <- cov.types[-j.strat]
# }
# j.trans <- grep("transformed",cov.types)
# #j.trans <- NULL
# j.cont <- grep("continuous",cov.types)
# cont.cov <- cov.names[j.cont]
# if (identical(covToTransform,"all"))
# covToTransform = cov.names[j.cont]
# tcov.names <- unique(c(cov.names[j.trans], setdiff(cont.cov,covToTransform)))
# cov.types <- gsub("transformed","",cov.types)
tcov.names <- NULL
covariates <- cov.info$covariate
cov.cat <- cov.names[cov.types == "categorical"]
covariates[cov.cat] <- lapply(covariates[cov.cat],as.factor)
# j.strat <- grep("stratification",cov.types)
# if (length(j.strat)>0)
# covariates <- covariates %>% select(-cov.names[j.strat])
#cov0 <- NULL
indvar <- mlx.getIndividualParameterModel()$variability$id
indvar[setdiff(param.names, paramToUse)] <- FALSE
cov.model <- mlx.getIndividualParameterModel()$covariateModel
r <- res <- r.cov0 <- list()
eps <- 1e-15
for (j in (1:n.param)) {
dj <- ind.dist[j]
nj <- names(dj)
if (indvar[j]) {
yj <- sp.df[nj]
if (tolower(dj) == "lognormal") {
yj <- log(yj+eps)
names(yj) <- paste0("log.",nj)
} else if (tolower(dj) == "logitnormal") {
yj <- log((yj+eps)/(1-yj+eps))
names(yj) <- paste0("logit.",nj)
} else if (tolower(dj) == "probitnormal") {
yj <- qnorm(yj)
names(yj) <- paste0("probit.",nj)
}
if (!is.null(eta))
etaj <- eta[paste0("eta_",nj)]
else
etaj <- NULL
if (length(covFix)>0) {
cmj <- cov.model[[nj]][covFix]
cov0 <- names(which(!cmj))
cov1 <- names(which(cmj))
} else {
cov0 <- cov1 <- NULL
}
# cov0 <- cov1 <- NULL
# browser()
# if (identical(weight,1))
# pwj <- 1
# else
pwj <- weight[nj,]
r[[j]] <- lm.all(nj, yj, etaj, covariates, tcov.names, pen.coef=pen.coef, nb.model=nb.model, pw=pwj, n.full=n.full,
direction=direction, steps=steps, p.max=p.max, cov0=cov0, cov1=cov1, iter=iter)
res[[j]] <- r[[j]]$res
r.cov0[[j]] <- r[[j]]$cov0
# res[[j]][,c("ll","df","criterion")] <- NULL
names(res[[j]]) <- gsub("log[.]","l",names(res[[j]]))
} else {
r[[j]] <- list(model="fixed")
res[[j]] <- "none"
}
r[[j]]$p.name <- nj
}
names(res) <- param.names
names(r.cov0) <- names(indvar)[which(indvar)]
e <- as.data.frame(lapply(r[indvar], function(x) {x$model$residuals}))
e.names <- unlist(lapply(r[indvar], function(x) {x$p.name}))
names(e) <- paste0("eta_",e.names)
if (!is.null(sp.df["id"]))
e <- cbind(sp.df["id"], e)
if (!is.null(sp.df["rep"]))
e <- cbind(sp.df["rep"], e)
if (length(correlation.model) >0 ) {
# en <- as.data.frame(scale(e))
for (ic in (1:length(correlation.model))) {
pic <- correlation.model[[ic]]
if (all(pic %in% e.names)) {
ipic <- match(pic, gsub("eta_","", names(e)))
epic <- e[,ipic]
gic <- solve(cov(epic))
jic <- match(pic, names(ind.dist))
# print(head(epic))
# if (1>2){
for (itk in 1:2) {
jk <- 0
for (j in jic) {
jk <- jk+1
dj <- ind.dist[j]
nj <- names(dj)
yj <- sp.df[nj]
pwj <- weight[nj,]
if (tolower(dj) == "lognormal") {
yj <- log(yj)
names(yj) <- paste0("log.",nj)
} else if (tolower(dj) == "logitnormal") {
yj <- log(yj/(1-yj))
names(yj) <- paste0("logit.",nj)
} else if (tolower(dj) == "probitnormal") {
yj[[1]] <- qnorm(yj[[1]])
names(yj) <- paste0("probit.",nj)
}
if (!is.null(eta))
etaj <- eta[paste0("eta_",nj)]
else
etaj <- NULL
if (length(covFix)>0) {
cmj <- cov.model[[nj]][covFix]
cov0 <- names(which(!cmj))
cov1 <- names(which(cmj))
} else {
cov0 <- cov1 <- NULL
}
ejc <- as.matrix(epic[,-jk])%*%matrix(gic[jk, -jk], ncol=1)/gic[jk, jk]
yjc <- yj + ejc
r[[j]] <- lm.all(nj, yjc, etaj, covariates, tcov.names, pen.coef=pen.coef, nb.model=nb.model, pw=pwj, n.full=n.full,
direction=direction, steps=steps, p.max=p.max, cov0=cov0, cov1=cov1, iter=iter)
res[[j]] <- r[[j]]$res
r.cov0[[j]] <- r[[j]]$cov0
# res[[j]][,c("ll","df","criterion")] <- NULL
#names(res[[j]]) <- gsub("log[.]","l",names(res[[j]]))
r[[j]]$p.name <- nj
e[paste0("eta_",nj)] <- epic[,jk] <- r[[j]]$model$residuals - ejc
}
# print(head(epic))
}
}
}
}
covariate.model <- mlx.getIndividualParameterModel()$covariateModel
covariate <- mlx.getCovariateInformation()$covariate
js <- 0
trs <- list()
tr0 <- NULL
for (k in (1:n.param)) {
if (!identical(res[[k]],"none")) {
# if (nrow(res[[k]]) < nb.model)
# res[[k]][(nrow(res[[k]])+1):nb.model, ] <- res[[k]][nrow(res[[k]]), ]
covariate.model[[k]][1:length(covariate.model[[k]])] <- FALSE
if (indvar[k]) {
ck <- attr(r[[k]]$model$terms,"term.labels")
if (length(ck)>0) {
for (j in (1:length(ck))) {
ckj <- ck[j]
if (identical(substr(ckj,1,4),"log.")) {
js <- js+1
ckj.name <- sub("log.","",ckj)
covkj <- covariate[[ckj.name]]
lckj <- paste0("l",ckj.name)
tr.str <- paste0(lckj,' = "log(',ckj.name,"/",signif(mean(covkj),digits=2),')"')
trs[[js]] <- paste0("lixoftConnectors::addContinuousTransformedCovariate(",tr.str,")")
tr0 <- unique(c(tr0,ckj.name))
#eval(parse(text=tr.str))
covariate.model[[k]][lckj] <- TRUE
} else {
covariate.model[[k]][ckj] <- TRUE
}
}
}
}
}
}
#res <- formatCovariateModel(getIndividualParameterModel()$covariateModel)
res <- formatCovariateModel(res)
return(list(model=covariate.model, residuals=e, res=res, add.covariate=trs,
sp=sp.df, tr0=tr0, r.cov0=r.cov0))
}
#-----------------------------------
lm.all <- function(ny, y, eta, x, tr.names=NULL, pen.coef=NULL, nb.model=NULL, pw=NULL, n.full=10,
direction='both', steps = 1000, p.max=1, cov0=NULL, cov1=NULL, iter=1) {
N <- length(unique(x$id))
nrep <- nrow(y)/N
if (p.max<=1) {
nx <- setdiff(names(x), c("id","rep"))
yc <- rowMeans(matrix(y[[1]], nrow=N))
# yc=colMeans(matrix(y[[1]], ncol=N))
xc <- x
# xc <- x[order(x$id),,drop=FALSE]
# xc <- x[seq(1,nrow(x),by=nrep),,drop=FALSE]
lm0 <- lm(yc ~1)
nxc <- setdiff(nx, cov0)
pjc <- NULL
for (nc in nxc) {
lmc <- lm(yc ~ xc[[nc]])
pc <- signif(anova(lm0, lmc)$`Pr(>F)`[2],4)
pjc <- c(pjc, pc)
# rc <- exp(pen.coef*(pw[nc]-1)/2)
# pjc <- c(pjc, pc*rc/(1 - pc + pc*rc))
}
pjc <- p.weight(pjc, pw[nxc], pen.coef)
names(pjc) <- nxc
if (!is.null(eta)) {
etac <- rowMeans(matrix(eta[[1]], nrow=N))
lm0 <- lm(etac ~1)
pjec <- NULL
for (nc in nxc) {
# if (length(etac) != length(xc[[nc]])) browser()
lmc <- lm(etac ~ xc[[nc]])
pc <- signif(anova(lm0, lmc)$`Pr(>F)`[2],4)
pjec <- c(pjec, pc)
# rc <- exp(pen.coef*(pw[nc]-1)/2)
# pjc <- c(pjc, pc*rc/(1 - pc + pc*rc))
}
pjec <- p.weight(pjec, pw[nxc], pen.coef)
names(pjec) <- nxc
# print(names(y))
# print(rbind(pjc, pjec))
pjc <- pmin(pjc, pjec, na.rm=T)
# if (names(y) %in% c("log.Cl", "log.V1")) pjc[1] <- 0
}
# if (length(which(pjc<p.max)) == 0)
# list.c <- which(pjc<max(pjc))
# else
# pjc <- p.adjust(pjc, method="BH")
pjc[names(which(mlx.getIndividualParameterModel()$covariateModel[[ny]]))] <- 0
list.c <- which(pjc>p.max)
# browser()
cov0 <- c(cov0, nxc[list.c])
direction <- ifelse(length(setdiff(nx, cov0))<=n.full,"full",direction)
# x <- x[, -which(names(x) %in% nx[list.c])]
} else
list.c <- NULL
x$id <- x$rep <- NULL
nx <- ncol(x)
l <- x
s <- rep("0:1", nx)
names(s) <- names(x)
j.num <- which(!sapply(x, is.factor))
j.num <- NULL
if (length(j.num)>0) {
if (length(j.num)==1)
j0 <- which(min(x[,j.num])>0)
else
j0 <- which(sapply(x[,j.num],min)>0)
j0.num <- j.num[j0]
s[j0.num] <- "0:2"
l[,j0.num] <- log(x[,j0.num])
names(l)[j0.num] <- paste0("log.",names(x)[j0.num])
# x[,j.num] <- scale(x[j.num], scale=TRUE)
# l[,j.num] <- scale(l[j.num], scale=TRUE)
l[,j.num] <- scale(l[j.num], scale=FALSE)
} else
j0.num <- vector(length=0)
#print(direction)
#############% stepAIC/stepBIC
if (direction != 'full') {
# create a data frame with all possible covariates (all transformations)
l_data = data.frame(y=y[[1]])
# not use non transformed covariates
j0.num = j0.num[!names(j0.num) %in% tr.names]
# not drop data frame to a vector!
l_data = cbind.data.frame(l_data,x,l[,j0.num,drop=FALSE])
# apply stepAIC
llk = tryCatch( {
# stepAIC forward, backward, both
f.sature <- "y ~ ."
if (length(cov0)>0)
f.sature <- paste0(f.sature, "-", paste(cov0,collapse="-"))
f.sature <- as.formula(f.sature)
model.sature=lm(f.sature,l_data)
f.cst <- "y ~ 1"
if (length(cov1)>0)
f.cst <- paste0(f.cst, "+", paste(cov1,collapse="+"))
f.cst <- as.formula(f.cst)
model.cst=lm(f.cst,l_data)
# model.sature=lm(y~.,l_data)
# model.cst=lm(y~1,data=l_data)
if(direction=='backward'){
lm.sat=mlx.stepAIC(model.sature, direction='backward', trace = FALSE, k=nrep*pen.coef,
scope=list(upper=model.sature, lower=model.cst), steps=steps, weight=pw)
}
else {
# lm.cst=mlx.stepAIC(model.cst, direction=direction, trace = FALSE, k=nrep*pen.coef,
# scope=list(upper=model.sature, lower=model.cst), steps=steps, weight=pw)
c.cur <- names(which(mlx.getIndividualParameterModel()$covariateModel[[ny]]))
f.cur <- "y ~ 1"
if (length(c.cur)>0)
f.cur <- paste0(f.cur, "+", paste(c.cur,collapse="+"))
f.cur <- as.formula(f.cur)
model.cur=lm(f.cur,l_data)
lm.cst=stepAIC(model.cur, direction=direction, trace=FALSE, k=nrep*pen.coef,
scope=list(upper=model.sature, lower=model.cst), steps=steps)
}
}
, error=function(e) {
print('Error in stepAIC')
return(-Inf) }
)
## after tryCatch
# Which covariates are used (written in 1-row data frame G)
Gnames = names(x)
G <- data.frame(matrix(0, ncol = length(Gnames), nrow = 1))
colnames(G) <- Gnames
# exact used covariates (with original names of cats), fill in G
usedcovariates = names(llk$model)[-1] #names(coef(llk)) 1, ]
G[1, usedcovariates] <- 1
# for (i in names(G)){
# if (i %in% usedcovariates ) G[1,i]=1 else G[1,i]=0
# if (paste0('log.',i) %in% usedcovariates ) G[1,i]=G[1,i]+2
# = 3 if both log and non-log are used
#}
ll = logLik(llk)/nrep
df = length(coef(llk))-1 # except Intercept
criterion = -2*ll + pen.coef*sum(pw[usedcovariates])
res <- data.frame(ll=round(ll,digits=3), df=df, criterion=round(criterion,digits=3))
j0.num <- j0.num[!(names(j0.num) %in% tr.names)]
if (length(j0.num)>0){
G[names(l)[j0.num]] <- 0
i2 <- (G[names(x)[j0.num]]==2)
G[names(l)[j0.num]][i2] <- 1
G[names(x)[j0.num]][i2] <- 0
}
res <- cbind(G==1, res)
# res <- cbind(G, res)
if (nb.model==1)
res[,c("ll","df","criterion")] <- NULL
else
res[2:nb.model,] <- res[1,]
row.names(res) <- 1:nrow(res)
# browser()
return(list(model=llk, res=res, cov0=cov0))
}
if (length(tr.names)>0){
j.newc <- which((names(x) %in% tr.names))
if (length(j.newc)>0)
s[j.newc] <- "0:1"
}
if (!is.null(cov0))
s[cov0] <- "0"
s <- paste(s,collapse=",")
s <- paste0("G <- expand.grid(", s, ")")
eval(parse(text=s))
if (length(tr.names)>0){
jc <- which(!(names(x) %in% tr.names))
c.names <- names(x)[jc]
for (k in 1:length(j.newc)) {
jk <- which(paste0("l",c.names)==tr.names[k])
if (length(jk)>0) {
j <- which(names(x) == c.names[jk])
tj <- which(names(x) == tr.names[k])
i0 <- which(G[,j]>0 & G[,tj]>0)
G <- G[-i0,]
}
}
}
names(G) <- names(x)
if (length(cov0)>0) {
i0 <- which(rowSums(G[cov0])==0)
G <- G[i0,]
}
if (length(cov1)>0) {
i1 <- which(rowSums(G[cov1]==1)==length(cov1))
G <- G[i1,]
}
ng <- nrow(G)
d <- ncol(G)
# print(c(ng,d))
ll <- df <- bic <- bic.cor <- NULL
corb <- log(iter^2/(iter^2+3))
#iop.mean <- ((iter ==2))
iop.mean <- F
if (iop.mean) {
yg <- colMeans(matrix(y[[1]], nrow=nrep))
xg <- x
if (nrow(l)>0)
lg <- l[seq(1,nrow(l),by=nrep),]
nrepg <- 1
# pen.coef <- pen.coef + log(nrep)
} else {
yg <- y[[1]]
xg <- x[rep(1:N,nrep),]
lg <- l
nrepg <- nrep
}
for (k in 1:ng) {
xk <- data.frame(y=yg)
Gk <- G[k,,drop=FALSE]
pwk <- pw[names(Gk)]
j1 <- which(Gk==1)
if (length(j1)>0)
xk[names(x)[j1]] <- xg[j1]
j2 <- which(Gk==2)
if (length(j2)>0)
xk[names(l)[j2]] <- lg[j2]
llk= tryCatch( {
lmk <- lm(y ~ ., data=xk)
logLik(lmk)[1]/nrepg }
, error=function(e) {
return(-Inf) }
)
dfk <- sum(Gk>0)
# dfk <- lmk$rank
# pen.covk <- 2*(sum(Gk*log((1-prior)/0.5)) + sum((1-Gk)*log(prior/0.5)))
# bick <- -2*llk + pen.coef*dfk + pen.covk
bick <- -2*llk + pen.coef*sum(Gk*pwk)
ll <- c(ll , llk)
df <- c(df, dfk)
bic <- c(bic , bick)
# bic.cor <- c(bic.cor , bick-corb*dfk)
bic.cor <- c(bic.cor , bick)
}
bic <- round(bic.cor, digits=3)
i0 <- rep(1,ng)
mG <- ncol(G)
for (k in seq_len(ng-1)) {
if (i0[k]==1) {
ik <- which(bic[(k):ng]==bic[k]) + k-1
sk <- .rowSums(G[ik,]==2, n=length(ik), m=mG)
ik0 <- ik[which(sk==0)]
if (length(ik0)==0)
ik0 <- ik[order(sk)[1]]
i0[ik] <- 0
i0[ik0] <- 1
i0[k] <- 1
}
}
res <- data.frame(ll=round(ll,digits=3), df=df, criterion=bic)
# print(res %>% group_by(df) %>% summarize(max(ll)))
# browser()
res <- res[i0==1,]
G <- G[i0==1,,drop=FALSE]
bic <- bic[i0==1]
#ill <- which(!duplicated(ll))
eval(parse(text=paste0(names(y)," <- y[[1]]")))
obic <- order(bic)
k.min <- obic[1]
# if (mG==1)
# Gkmin <- G[k.min]
# else
Gkmin <- G[k.min,]
j1 <- which(Gkmin==1)
j2 <- which(Gkmin==2)
if (length(j1)>0) {
for (k in (1:length(j1)))
eval(parse(text=paste0(names(x)[j1[k]]," <- rep(x[[j1[k]]], nrep)")))
}
if (length(j2)>0) {
for (k in (1:length(j2)))
eval(parse(text=paste0(names(l)[j2[k]]," <- l[[j2[k]]]")))
}
list.x <- c("1",names(x)[j1],names(l)[j2])
form1 <- paste0(names(y), "~", paste(list.x, collapse = "+"))
eval(parse(text=paste0("lm.min <- lm(",form1,")")))
lm.min$covsel=Gkmin
nb.model0 <- min(nb.model, length(bic))
res <- res[obic[1:nb.model0],]
G <- G[obic[1:nb.model],,drop=FALSE]
j0.num <- j0.num[!(names(j0.num) %in% tr.names)]
if (length(j0.num)>0){
G[names(l)[j0.num]] <- 0
i2 <- (G[names(x)[j0.num]]==2)
G[names(l)[j0.num]][i2] <- 1
G[names(x)[j0.num]][i2] <- 0
}
if (nb.model>nb.model0)
res[(nb.model0+1):nb.model,c("ll","df","criterion")] <- NA
res <- cbind(G==1, res)
if (nb.model==1)
res[,c("ll","df","criterion")] <- NULL
row.names(res) <- 1:nrow(res)
return(list(model=lm.min, res=res, cov0=cov0))
}
#######################################################
mlx.stepAIC<-function (object, scope, scale = 0, direction = c("both", "backward","forward"),
trace = 1, keep = NULL, steps = 1000, use.start = FALSE, k = 2, weight=NULL)
{
mydeviance <- function(x) {
dev <- deviance(x)
if (!is.null(dev))
dev
else extractAIC(x, k = 0)[2L]
}
cut.string <- function(string) {
if (length(string) > 1L)
string[-1L] <- paste("\n", string[-1L], sep = "")
string
}
re.arrange <- function(keep) {
namr <- names(k1 <- keep[[1L]])
namc <- names(keep)
nc <- length(keep)
nr <- length(k1)
array(unlist(keep, recursive = FALSE), c(nr, nc), list(namr, namc))
}
step.results <- function(models, fit, object, usingCp = FALSE) {
change <- sapply(models, "[[", "change")
rd <- sapply(models, "[[", "deviance")
dd <- c(NA, abs(diff(rd)))
rdf <- sapply(models, "[[", "df.resid")
ddf <- c(NA, abs(diff(rdf)))
AIC <- sapply(models, "[[", "AIC")
heading <- c("Stepwise Model Path \nAnalysis of Deviance Table",
"\nInitial Model:", deparse(formula(object)), "\nFinal Model:",
deparse(formula(fit)), "\n")
aod <- if (usingCp)
data.frame(Step = change, Df = ddf, Deviance = dd,
`Resid. Df` = rdf, `Resid. Dev` = rd, Cp = AIC,
check.names = FALSE)
else data.frame(Step = change, Df = ddf, Deviance = dd,
`Resid. Df` = rdf, `Resid. Dev` = rd, AIC = AIC,
check.names = FALSE)
attr(aod, "heading") <- heading
class(aod) <- c("Anova", "data.frame")
fit$anova <- aod
fit
}
extractAIC.getpenalisationweight<-function(object, weight, k){
variableyn<-weight
variableyn[]<-0
variableyn[attr(terms(object),"term.labels")]<-1
penalweight<- sum(variableyn*(weight-1)*k)
penalweight
}
extractAIC.updateAICweight<-function(object,objaod,weight,k,drop=TRUE){
objaod[,"AIC"]<-objaod[,"AIC"] + extractAIC.getpenalisationweight(object,weight,k)
no <- rownames(objaod)[2:dim(objaod)[1]]
if (drop)
objaod[no,"AIC"] <- objaod[no,"AIC"] - (weight[no]-1)*k
else
objaod[no,"AIC"] <- objaod[no,"AIC"] + (weight[no]-1)*k
objaod
}
Terms <- terms(object)
object$formula <- Terms
if (inherits(object, "lme"))
object$call$fixed <- Terms
else if (inherits(object, "gls"))
object$call$model <- Terms
else object$call$formula <- Terms
if (use.start)
warning("'use.start' cannot be used with R's version of 'glm'")
md <- missing(direction)
direction <- match.arg(direction)
backward <- direction == "both" | direction == "backward"
forward <- direction == "both" | direction == "forward"
if (missing(scope)) {
fdrop <- numeric()
fadd <- attr(Terms, "factors")
if (md)
forward <- FALSE
}
else {
if (is.list(scope)) {
fdrop <- if (!is.null(fdrop <- scope$lower))
attr(terms(update.formula(object, fdrop)), "factors")
else numeric()
fadd <- if (!is.null(fadd <- scope$upper))
attr(terms(update.formula(object, fadd)), "factors")
}
else {
fadd <- if (!is.null(fadd <- scope))
attr(terms(update.formula(object, scope)), "factors")
fdrop <- numeric()
}
}
models <- vector("list", steps)
if (!is.null(keep))
keep.list <- vector("list", steps)
n <- nobs(object, use.fallback = TRUE)
fit <- object
bAIC <- extractAIC(fit, scale, k = k) #Was +...
edf <- bAIC[1L]
bAIC <- bAIC[2L]+extractAIC.getpenalisationweight(fit,weight,k)
if (is.na(bAIC))
stop("AIC is not defined for this model, so 'stepAIC' cannot proceed")
if (bAIC == -Inf)
stop("AIC is -infinity for this model, so 'stepAIC' cannot proceed")
nm <- 1
Terms <- terms(fit)
if (trace) {
cat("Start: AIC=", format(round(bAIC, 2)), "\n", cut.string(deparse(formula(fit))),
"\n\n", sep = "")
utils::flush.console()
}
models[[nm]] <- list(deviance = mydeviance(fit), df.resid = n -
edf, change = "", AIC = bAIC)
if (!is.null(keep))
keep.list[[nm]] <- keep(fit, bAIC)
usingCp <- FALSE
while (steps > 0) {
steps <- steps - 1
AIC <- bAIC
ffac <- attr(Terms, "factors")
if (!is.null(sp <- attr(Terms, "specials")) && !is.null(st <- sp$strata))
ffac <- ffac[-st, ]
scope <- factor.scope(ffac, list(add = fadd, drop = fdrop))
aod <- NULL
change <- NULL
if (backward && length(scope$drop)) {
aod <- dropterm(fit, scope$drop, scale = scale, trace = max(0, trace - 1), k = k)#Was +...
aod<-extractAIC.updateAICweight(fit,aod,weight,k,drop=TRUE)
rn <- row.names(aod)
row.names(aod) <- c(rn[1L], paste("-", rn[-1L], sep = " "))
if (any(aod$Df == 0, na.rm = TRUE)) {
zdf <- aod$Df == 0 & !is.na(aod$Df)
nc <- match(c("Cp", "AIC"), names(aod))
nc <- nc[!is.na(nc)][1L]
ch <- abs(aod[zdf, nc] - aod[1, nc]) > 0.01
if (any(is.finite(ch) & ch)) {
warning("0 df terms are changing AIC")
zdf <- zdf[!ch]
}
if (length(zdf) > 0L)
change <- rev(rownames(aod)[zdf])[1L]
}
}
if (is.null(change)) {
if (forward && length(scope$add)) {
aodf <- addterm(fit, scope$add, scale = scale, trace = max(0, trace - 1), k = k)#Was +...
aodf<-extractAIC.updateAICweight(fit, aodf, weight, k, drop=FALSE)
rn <- row.names(aodf)
row.names(aodf) <- c(rn[1L], paste("+", rn[-1L],
sep = " "))
aod <- if (is.null(aod))
aodf
else rbind(aod, aodf[-1, , drop = FALSE])
}
attr(aod, "heading") <- NULL
if (is.null(aod) || ncol(aod) == 0)
break
nzdf <- if (!is.null(aod$Df))
aod$Df != 0 | is.na(aod$Df)
aod <- aod[nzdf, ]
if (is.null(aod) || ncol(aod) == 0)
break
nc <- match(c("Cp", "AIC"), names(aod))
nc <- nc[!is.na(nc)][1L]
o <- order(aod[, nc])
if (trace) {
print(aod[o, ])
utils::flush.console()
}
if (o[1L] == 1)
break
change <- rownames(aod)[o[1L]]
}
usingCp <- match("Cp", names(aod), 0) > 0
fit <- update(fit, paste("~ .", change), evaluate = FALSE)
fit <- eval.parent(fit)
nnew <- nobs(fit, use.fallback = TRUE)
if (all(is.finite(c(n, nnew))) && nnew != n)
stop("number of rows in use has changed: remove missing values?")
Terms <- terms(fit)
bAIC <- extractAIC(fit, scale, k = k)#Was +...
edf <- bAIC[1L]
bAIC <- bAIC[2L]+extractAIC.getpenalisationweight(fit,weight,k)
if (trace) {
cat("\nStep: AIC=", format(round(bAIC, 2)), "\n",
cut.string(deparse(formula(fit))), "\n\n", sep = "")
utils::flush.console()
}
if (bAIC >= AIC + 1e-07)
break
nm <- nm + 1
models[[nm]] <- list(deviance = mydeviance(fit), df.resid = n -
edf, change = change, AIC = bAIC)
if (!is.null(keep))
keep.list[[nm]] <- keep(fit, bAIC)
}
if (!is.null(keep))
fit$keep <- re.arrange(keep.list[seq(nm)])
step.results(models = models[seq(nm)], fit, object, usingCp)
}
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.