Nothing
frailty1.vs <- function(formula,model,penalty,data,B="null",v="null",alpha="null",tun1="null",tun2="null",varfixed=FALSE,varinit=0.1){
Call <- match.call()
mc <- match.call()
indx <- match(c("formula", "data"),
names(Call), nomatch = 0)
if (indx[1] == 0)
stop("A formula argument is required")
temp <- Call[c(1, indx)]
temp[[1]] <- as.name("model.frame")
special <- c("strata", "cluster")
temp$formula <- terms(subbars(formula), special)
m <- eval(temp)
Terms <- attr(m, "terms")
Y <- model.extract(m, "response")
temp <- Call[c(1, indx)]
temp[[1]] <- as.name("model.frame")
special <- c("strata", "cluster")
temp$formula <- terms(formula, special)
Terms <- temp[[2]]
formula1 <- paste(paste(Terms[[2]][[2]], Terms[[3]], sep = "~")[[2]],
paste(Terms[[3]])[3], sep = "+")
formula1 <- formula(formula1)
fr <- FrailtyFrames(mc, formula1, contrasts)
namesX <- names(fr$fixef)
namesX <- namesX[-1]
namesY <- names(fr$mf)[1]
FL <- HGLMFactorList(formula1, fr, 0L, 0L)
namesRE <- FL$namesRE
y <- matrix(Y[, 1], length(fr$Y), 1)
x <- fr$X
z <- FL$Design
n <- nrow(x)
p <- ncol(x)
x1 <- x[1:n, 2:p]
x2 <- matrix(x1, n, p - 1)
x <- x2
n <- nrow(x)
p <- ncol(x)
nrand <- length(z)
q <- rep(0, nrand)
for (i in 1:nrand) q[i] <- dim(z[[i]])[2]
del <- matrix(0, n, 1)
del[, 1] <- censor <- Y[, 2]
SS <- FL$Subject
res1 <- FrailtyMakeData(y, x, del, z)
y <- as.numeric(res1[1][[1]])
x <- as.matrix(res1[2][[1]])
del <- as.numeric(res1[3][[1]])
z <- res1[4][[1]]
Mi <- as.matrix(res1[5][[1]])
idx2 <- res1[6][[1]]
t2 <- as.numeric(res1[7][[1]])
di <- as.numeric(res1[8][[1]])
if(is.null(v)){
v<-NULL
for (i in 1:nrand) v[[i]]<-rep(0, q[i])
}
if(is.null(alpha)){
alpha<-rep(0.1, nrand)
#alpha<-c(0.026, 1.036)
}
if (varfixed==TRUE){alpha<-rep(max(varinit, 1e-06),nrand)}
if(length(B)==1){
if(is.null(B)){
B<-as.numeric(frailtyHL(formula,data=data)$FixCoef[,1])
}
}
#print(z)
if((length(tun1)!=1)|(length(tun2)!=1)){
if (model=="lognorm"){
result<-lognorm1.BIC(X=x,B,Z=z,v,penalty,alpha,tun1,tun2,di,Mi,idx2,del,varfixed)
}
#if (model=="gamma"){
#result<-gammapv.BIC(X=x,B,Z=z,v,penalty,alpha,tun1,tun2,di,Mi,idx2,del)
#}
}
else{
if (model=="lognorm"){
result<-lognorm1.iter(X=x,B,Z=z,v,penalty,alpha,tun1,tun2,di,Mi,idx2,del,varfixed)
}
#if (model=="gamma"){
#result<-gammapv.iter(X=x,B,Z=z,v,penalty,alpha,tun1,tun2,di,Mi,idx2,del)
#}
}
print("Result of variable selection in frailty model")
print("==Fitted model==")
print(paste("model :", model))
print(paste("penalty :", penalty))
print("formula :")
print(formula)
cn<-c("did not converge")
if(result$convergence==1){
cn<-c("converge")
}
print(cn)
print("==Fixed coefficients==")
print_mat<-cbind(result$beta,result$se)
rownames(print_mat)<-namesX
colnames(print_mat)<-c("Estimate","Std. Error")
print(print_mat)
print("==Dispersion parameter==")
print(result$alpha)
print("==Tuning parameter==")
print(result$tuning_parameter)
print("==BIC==")
print(as.numeric(result$BIC))
result$res1<-res1
return(result)
}
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.