# R/concreg.r In concreg: Concordance Regression

#### Documented in cindexcoef.concregconcregconfint.concregvcov.concreg

```useFP <- function
(
obj,
data
)
### alternative idea to get plot fp terms etc ... missing a lot ...
### 2008-11
{
browser()
obj.full <- decomposeSurv(obj\$formula, data, sort=FALSE)
## -> \$fac,resp,mm1,NTDE,timedata,timeind,NFP,fpnames,fpind,PTcoefs,covnames,ind,offset.values

## take obj.full\$mm1[, obj\$ind]

CODE <- paste("(I(PT(\\1)) + powM2(PT(\\1)) + powM1(PT(\\1)) + powM0.5(PT(\\1)) + log(PT(\\1)) + ",
"sqrt(PT(\\1)) + pow2(PT(\\1)) + pow3(PT(\\1)) + ",
"RI(PT(\\1)) + RpowM2(PT(\\1)) + RpowM1(PT(\\1)) + RpowM0.5(PT(\\1)) + ",
"Rlog(PT(\\1)) + Rsqrt(PT(\\1)) + Rpow2(PT(\\1)) + Rpow3(PT(\\1)) )")
sub3 <- gsub("fp\\(([a-zA-Z0-9]*)\\)", CODE, obj\$formula[3])
formula <- as.formula(paste(as.character(obj\$formula)[2], "~", sub3))

## define simple transformations
powM2 <- function(z) z^(-2)
powM1 <- function(z) z^(-1)
powM0.5 <- function(z) z^(-0.5)
pow2 <- function(z) z^2
pow3 <- function(z) z^3

## define repeated powers
RI <- function(z) z * log(z)
RpowM2 <- function(z) z^(-2) * log(z)
RpowM1 <- function(z) z^(-1) * log(z)
RpowM0.5 <- function(z) z^(-0.5) * log(z)
Rlog <- function(z) log(z) * log(z)
Rsqrt <- function(z) sqrt(z) * log(z)
Rpow2 <- function(z) z^2 * log(z)
Rpow3 <- function(z) z^3 * log(z)

##
browser()

}

fp.scale <- function
(
x
)
### taken from package <mfp>
### 2008-06
{
scale <- 1
shift <- 0
if (min(x) <= 0) {
z <- diff(sort(x))
shift <- min(z[z > 0]) - min(x)
shift <- ceiling(shift * 10)/10
}
range <- mean(x + shift)
scale <- 10^(sign(log10(range)) * round(abs(log10(range))))

list(
shift=shift,
scale=scale
)
}

decomposeSurv <- function
(
formula,
data,
sort=FALSE,
offset=NULL
)
### decomposes complex survival formula
### trans: I(), powM2, powM1, powM0.5, sqrt, pow2, pow3, log,
### 2008-04
{
orig.formula <- formula

## expand formula if needed:
repeat {
terms <- terms(formula, "fp", data=data)
fac <- attr(terms, "factors")

needed <- rownames(fac)[!rownames(fac) %in% colnames(fac)][-1]

if(length(needed) == 0) break

formula <- as.formula(paste(as.character(formula)[2], "~",
as.character(formula)[3], "+",
paste(needed, sep="+")))
}

## determine position of <fp> terms (of all orders) in
## not expanded formula
spec <- untangle.specials(terms, "fp", order=1:3)
fpnames <- gsub("fp\\(([a-zA-Z0-9]*)\\)", "\\1", spec\$vars) # c("a", "b")
fpiden <- paste("(", fpnames, ")", sep="") # c("(a)", "(b)")

formulaOrig <- formula
## replace fp() by the 8x2 possible terms (in formula[3])
## attention: fp-variables must consist of a-z,A-Z,0-9, otherwise
## must be made modifications here (gsub) ...
## FURTHER: 1st term is the linear one
## FURTHER: 1st half (=8) are the usual terms, 2nd half (9-16) are the repeated powers
## FURTHER: problems possible if one variable name is part of the other, e.g. ABC & ABCDE
##CODE <- paste("(I(\\1) + powM2(\\1) + powM1(\\1) + powM0.5(\\1) + log(\\1) + sqrt(\\1) + pow2(\\1) + pow3(\\1) + ",
##              "  RI(\\1) + RpowM2(\\1) + RpowM1(\\1) + RpowM0.5(\\1) + Rlog(\\1) + Rsqrt(\\1) + Rpow2(\\1) + Rpow3(\\1) )")
CODE <- paste("(I(PT(\\1)) + powM2(PT(\\1)) + powM1(PT(\\1)) + powM0.5(PT(\\1)) + log(PT(\\1)) + ",
"sqrt(PT(\\1)) + pow2(PT(\\1)) + pow3(PT(\\1)) + ",
"RI(PT(\\1)) + RpowM2(PT(\\1)) + RpowM1(PT(\\1)) + RpowM0.5(PT(\\1)) + ",
"Rlog(PT(\\1)) + Rsqrt(PT(\\1)) + Rpow2(PT(\\1)) + Rpow3(PT(\\1)) )")
sub3 <- gsub("fp\\(([a-zA-Z0-9]*)\\)", CODE, formulaOrig[3])
formula <- as.formula(paste(as.character(formula)[2], "~", sub3))

## define simple transformations
powM2 <- function(z) z^(-2)
powM1 <- function(z) z^(-1)
powM0.5 <- function(z) z^(-0.5)
pow2 <- function(z) z^2
pow3 <- function(z) z^3

## define repeated powers
RI <- function(z) z * log(z)
RpowM2 <- function(z) z^(-2) * log(z)
RpowM1 <- function(z) z^(-1) * log(z)
RpowM0.5 <- function(z) z^(-0.5) * log(z)
Rlog <- function(z) log(z) * log(z)
Rsqrt <- function(z) sqrt(z) * log(z)
Rpow2 <- function(z) z^2 * log(z)
Rpow3 <- function(z) z^3 * log(z)

## pretransformation function
PT <- function(z) {
obj <- fp.scale(z)
(z + obj\$shift) / obj\$scale }

## construct 3-col response:
resp <- model.extract(model.frame(formula, data = data), "response")
if(is.null(dim(resp))) 	 resp<-cbind(resp-min(min(resp-0.01),0), rep(1,length(resp)))
if(ncol(resp) == 2)      resp <- cbind(start=rep(0, nrow(resp)), resp)

## sortieren nach STOPzeit und -Cens
if(sort) {
sort <- order(resp[, 2],  -resp[, 3])
data <- data[sort, , drop=FALSE]
resp <- resp[sort, ]
}

mm <- model.matrix(formula, data = data) ## Model-Matrix
mm1 <- mm[, -1, drop=FALSE]	# w/o intercept

## *** RETREIVE PRETRANS COEFS
sub3 <- gsub("fp\\(([a-zA-Z0-9]*)\\)", "PT(\\1)", formulaOrig[3])
formulaPT <- as.formula(paste(as.character(formulaOrig)[2], "~", sub3))
## varied pretransformation function to get coefficients:
PT <- function(z) {
obj <- fp.scale(z)
c(shift=obj\$shift, scale=obj\$scale, rep(-999, length(z) - 2)) }
mmPT <- model.matrix(formulaPT, data = data)
PTcoefs <- mmPT[1:2, !is.na(mmPT[1,]) & (mmPT[3,]==-999), drop=FALSE]
rownames(PTcoefs) <- c("shift", "scale")
## *** END

## offset ..
if(length(offset) != 0)
offset.values <- offset
else
offset.values <- NA

terms <- terms(formula, "fp", data=data)
fac <- attr(terms, "factors")
labels <- attr(terms, "term.labels")

## splits by special chars
f <- function(str)
for(chars in c("(", ")", ":", " ", ",", "*", "^"))
str <- unlist(strsplit(str, split=chars, fixed=TRUE))

rowSplit <- sapply(rownames(fac), f, simplify=FALSE)	# splitted effects
stopName <- tail(rowSplit[[1]], 2)[1]	# name of stoptime
rowInter <- unlist(lapply(rowSplit[-1], function(z) any(z == stopName)))
##  rowOffset <- unlist(lapply(rowSplit[-1], function(z) any(z == "offset")))
##  rowOffset

fac <- fac[-1, , drop=FALSE]	# omit Surv

colSplit <- lapply(colnames(fac), f)
colInter <- unlist(lapply(colSplit, function(z) any(z == stopName)))

nTimes <- colSums(fac[rowInter, , drop=FALSE])
nFac   <- colSums(fac[!rowInter, , drop=FALSE])

inters <- (nFac>0) & (nTimes>0)
NTDE <- sum(inters)

timedata <- matrix(0, nrow(data), 0)
timeind <- c()

## loop for (time x effect)
for(i in which(inters)) {
## search pure time:
ind <- (colSums(fac[rowInter, i] != fac[rowInter, , drop=FALSE]) == 0) & (nFac==0)
timedata <- cbind(timedata, mm1[, ind, drop=FALSE])

## search pure effect:
ind <- (colSums(fac[!rowInter, i] != fac[!rowInter, , drop=FALSE]) == 0) & (nTimes == 0)
timeind <- c(timeind, which(ind[!colInter]))
}
mm1 <- mm1[, !colInter, drop=FALSE]

covnames <- c(colnames(mm1),
paste(colnames(timedata), colnames(mm1)[timeind], sep=":")
)
alter <- c(colnames(mm1),
paste(colnames(mm1)[timeind], colnames(timedata), sep=":")
)

## indicator to identify the original formula:
ind <- covnames %in% colnames(attr(terms(orig.formula, "fp", data=data), "factors")) |
alter %in% colnames(attr(terms(orig.formula, "fp", data=data), "factors"))

## FP indicator matrix (1:16)
NFP <- length(fpnames)
fpind <- matrix(0, NFP, length(covnames), dimnames=
list(fpnames, covnames))
for(i in seq(length=NFP)) {
inds <- grep(fpiden[i], covnames)
## when interactions occur, 1:16,1:16 instead of 1:32 is needed
fpind[i, inds] <- ((seq(along=inds) - 1) %% 16) + 1
}

## return object
list(
fac=fac,                    # factor matrix ..
resp=resp,                  # N x 3 - response matrix
mm1=mm1,                    # model matrix without time effects

NTDE=NTDE,                  # number time dep. effects
timedata=timedata,          # matrix with time functions as columns
timeind=timeind, 		       # indicator of time-dependend effects

NFP=NFP,                    # number of frac.polys
fpnames=fpnames,            # names of the variables used as fractional poly's
fpind=fpind,                # matrix with frac.polyn. in each row, numbered by 1 to 5
PTcoefs=PTcoefs,            # coefficients of pretransformations

covnames=covnames,          # names of covariates
ind=ind,			               # indicator vector:
# which terms are really part of the formula
offset.values=offset.values # offset values
)
}

concreg <- function
(
formula=attr(data, "formula"),         # formula
data=sys.parent(),                     #
id=NULL,                               # identifier: numeric or character or factor
normalize=TRUE,
scale.weights=1,
offset=NULL,
alpha=0.05,                            # confidence limit
maxit=50,                              # max. iterations
maxhs=5,                               # half steps
epsilon=1e-6,                          #
maxstep=2.5,                           #
x=TRUE,                                # for output
y=TRUE,                                # for output
print=TRUE,                            # print fitting information on screen
c.risk=NULL,                           # Vector encoding competing risk: 0=censored, 1=event, 2=competing risk
strata.var=NULL,                       # strata variable name
trunc.weights=1,                       # quantile for weight truncation: all weights greater than that quantile will be truncated to that value
npar=FALSE,                            # nonparametric estimation?
# breslow=NA, # righthand formula, if breslow weighted terms, e.g. ~ A + A:C
# prentice=NA, # righthand formula, if prentice weighted terms, e.g. ~ A + A:C
# taroneware=NA, # righthand formula, if tarone-ware weighted terms, e.g. ~ A + A:C
# censcorr=FALSE,
# robust=FALSE,
# jack=FALSE,
...                                    # N -> breslow; km -> prentice
)
### by MP und GH, 2008
### Surv-Objekt entweder usual  (z.B. Surv(time,event)~A+B+C+D)
### event: 2=competing risk, 1=dead/event, 0=censored
{
## evt. use of abbrev. parameter names
#        l <- list(...)
#        if("N" %in% names(l)) breslow <- l\$N
#        if("km" %in% names(l)) prentice <- l\$km
alpha.fp=c(0.20, 0.05)
fp.iter=10                            # maximum number of iterations of large <fp> loop
n <- nrow(data)
#        pl <- FALSE
#        robcov <- if(robust) 1 else 0

## generate or reorder id's such that values are within 1:n
if (is.null(id)) { id <- 1:n } else id <- as.numeric(as.factor(id))
maxid <- max(id)

## here only ONCE the full model matrix is spanned with all possible fp-effects
obj.full <- decomposeSurv(formula, data, sort=FALSE, offset)

# change Daniela Competing Risk
if (is.null(c.risk)) {
crisk <- obj.full\$resp[,3]
obj.full\$resp <- cbind(obj.full\$resp, crisk)

}

if (is.null(c.risk)==FALSE) {
crisk <- obj.full\$resp[,3]
c.risk[c.risk==1] <- 2
crisk[crisk==0] <- c.risk[crisk==0]
obj.full\$resp <- cbind(obj.full\$resp, crisk)
}

## stratagruppen bilden
if (is.null(strata.var)[1]) {
obj.full\$stratum <- rep(1, n)
}

if (is.null(strata.var)[1] == FALSE) {
mult.strata <- 1
if (length(strata.var) > 1)  {
for (l in 1:length(strata.var)) { mult.strata <- c(mult.strata, mult.strata[length(mult.strata)]*10) }
}

daten.tmp <- data.frame(data[,strata.var])
for (l in 1:length(strata.var)) {
anz.strata <- length(table(data[,strata.var[l]]))

for (l1 in 1:anz.strata)  {                                               # Strata neu durchnummerieren
daten.tmp[data[,strata.var[l]] == as.numeric(names(table(data[,strata.var[l]])))[l1], l] <- l1
}

daten.tmp[, l] <- daten.tmp[, l]*mult.strata[l]
}

if (length(strata.var)==1)  { daten.tmp\$tmp <- daten.tmp  }
if (length(strata.var)> 1)  { daten.tmp\$tmp <- apply(daten.tmp[, strata.var], 1, sum) }
tab <- table(daten.tmp\$tmp)

#    if (sum(tab==1)!=0) { stop("number of observations is in at least 1 stratum 1")}     #! ev. einfach diese Daten weglassen
anz.strata <- length(tab)
for (l2 in 1:anz.strata)  {
daten.tmp\$stratum[daten.tmp\$tmp == as.numeric(names(tab)[l2])] <- l2
}
obj.full\$stratum <- daten.tmp\$stratum
}

max.strata <- max(obj.full\$stratum)

# change Daniela Competing Risk Weights (getrennt nach Strata)
# geht nicht so einfach, Problem mit ties (gleiche Zeiten nur 1 mal in der Liste)
if(!is.null(c.risk))   #### Georg 110517
{
vgl <- G <- data.frame(cbind(obj.full\$resp[,2], obj.full\$stratum, -99))
G <- unique(G)
for (i.strata in 1: max.strata) {
g.tmp <- rep(NA, times=sum(obj.full\$stratum==i.strata))
g.tmp[obj.full\$resp[obj.full\$stratum==i.strata,3]==0] <- 1
g.tmp[obj.full\$resp[obj.full\$stratum==i.strata,3]==1 | obj.full\$resp[obj.full\$stratum==i.strata,3]==2] <- 0

G[G\$X2==i.strata,3] <- survfit(Surv(obj.full\$resp[obj.full\$stratum==i.strata,2], g.tmp)~1)\$surv
}

if (nrow(G) != nrow(vgl))  {                        # Bei Bindungen innerhalb der Strata
G.n <- 1
for (i.n in 1:n)  {
if (sum(vgl[i.n, 1:2] == G[G.n, 1:2])==2) {
vgl[i.n, 3] <- G[G.n,3]
G.n <- G.n+1
}
}

vgl[order(vgl[,3])[1:sum(vgl[,3] == -99)],3] <- vgl[order(vgl[,3])[1:sum(vgl[,3] == -99)]-1,3]
G <- vgl
}
G <- G[,3]

} else G<-rep(1,nrow(obj.full\$resp))     #### Georg 110517
# new time variable for competing risk
time.crisk <- obj.full\$resp[,2]
time.crisk[obj.full\$resp[,"crisk"]==2] <- max(obj.full\$resp[,2])+1
obj.full\$resp <- cbind(obj.full\$resp, time.crisk)

# calculate weights
W <- concreg.wei(resp=obj.full\$resp, max.strata=max.strata, stratum=obj.full\$stratum, trunc.weights=trunc.weights)

obj <- obj.full

kk <- ncol(obj.full\$mm1) # should be k2 - NTDE
obj\$mm1 <- obj.full\$mm1[, obj\$ind[1:kk], drop=FALSE]
obj\$covnames <- obj.full\$covnames[obj\$ind]

obj\$timeind  <- obj.full\$timeind[obj\$ind[-(1:kk)]]                      #! weg?
obj\$timedata <- obj.full\$timedata[, obj\$ind[-(1:kk)], drop=FALSE]       #! weg?
## re-index \$timeind
obj\$timeind <- match(obj\$timeind, (1:kk)[obj\$ind[1:kk]])                #! weg?
NTDE <- obj\$NTDE <- length(obj\$timeind)                                 #! weg?

ind.offset <- sum(length(offset) != 0)

k <- ncol(obj\$mm1)    # number covariates w/o time-dep effects          # Anzahl Variablen
k2 <- k + NTDE                                                          #! weg?
if (npar) nonpar<-1
else nonpar<-0

PARMS <- c(n, k, 0, maxit, maxhs, maxstep, epsilon,  0, 0, 0, 0, 0, 0, max.strata, maxid, ind.offset, nonpar)

# alles sortieren (obj, W, G, id ist sortiert)
ord <- order(obj.full\$stratum, obj.full\$resp[,"time.crisk"],(1-obj.full\$resp[,3]))
obj\$mm1     <- obj\$mm1[ord,]
obj\$resp    <- obj\$resp[ord,]
obj\$stratum <- obj\$stratum[ord]
G           <- G[ord]
W           <- W[ord]
id          <- id[ord]

##   if (offset) {
##    IOARRAY[1,1]<-0    # first variable is offset
##    IOARRAY[2,1]<-Z.sd[1]    # first variable is offset
##   }
## **************** fit model *********************
value0 <- concreg.fit(obj=obj, id=id, W=W, G=G, PARMS=PARMS, npar=npar)            #! Aufruf der Funktion mit dem Fortran Aufruf

#        cov.ls <- value0\$cov.ls
#        cov.lw <- cov.j <- dfbeta.resid <- NULL

#        if(robust) {
#                cov.lw <- matrix(value0\$outtab[(k2+4):(3+2*k2), ], ncol=k2) / value0\$ZxZ
#                dfbeta.resid <- value0\$dfbetaresid / matrix(value0\$Z.sd, maxid, k2, byrow=TRUE)
#                covs <- cov.lw
#                cov.method <- "Lin-Wei"
#
#        } else if(jack) {
#                dfbeta.resid <- matrix(0, maxid, k2)
#                nc <- ncol(value0\$cards)
#                for(iid in 1:maxid) {
#                        cardsJ <- value0\$cards[id != iid,]
#                        cardsJ[id[id!=iid]>iid,nc] <- cardsJ[id[id!=iid]>iid,nc]-1
#                        n.jack <- nrow(cardsJ)
#                        PARMS <- c(n.jack, k, 0, maxit, maxhs,maxstep,epsilon, 0, 0, 0, 0, 0, W\$NGV, NTDE, maxid-1, ind.offset)
#                        valueJ <- coxphw.fit(obj, id, W\$weights, PARMS, cardsJ) # fit model jack(i) *******
#                        dfbeta.resid[iid,] <- value0\$coef.orig - valueJ\$coef.orig
#                }
#                cov.j <- (maxid-1) / maxid * crossprod(dfbeta.resid) / value0\$ZxZ
#                dfbeta.resid <- dfbeta.resid / matrix(value0\$Z.sd, maxid, k2, byrow=TRUE)
#                covs <- cov.j
#                cov.method <- "Jackknife"
#
#        } else {
#                cov.method <- "Lin-Sasieni"
#                covs <- cov.ls
#        }

vars <- as.matrix(diag(value0\$cov.rob))
#        dimnames(vars) <- list(obj\$covnames, obj\$covnames)

## DECIDE THE MODEL: USE PVALS ##############
probs <- 1 - pchisq((value0\$coefs^2/vars), 1)

## ########## NOW ONLY FINAL MODEL IS CONSIDERED ############
names(value0\$coefs) <- obj\$covnames
if(value0\$outpar[10]>=maxit)
cat("No convergence attained in ", value0\$outpar[10], " iterations.\n", sep="")

Means <- colMeans(value0\$mmm)

## return object
fit <- list(coefficients = value0\$coefs,     # coefficients of the fit
cards    = value0\$cards,         #
parms    = value0\$outpar,
ioarray  = value0\$outtab,
#                    dfbeta.resid = dfbeta.resid,
alpha    = alpha,                # significance level
var      = value0\$cov.rob,                 # covariance matrix
df       = k2,                   # degrees of freedom (k + NTDE)
iter     = value0\$outpar[10],
method.ties = "no",         #
n = n,                           # number observations
##terms = terms(formula),
y = obj\$resp,                    # responses
formula = formula,               # original formula
exit.code=value0\$outpar[8],

#                  fpind=obj\$fpind,
#                  PTcoefs=obj.full\$PTcoefs,                                   #! weg?
#                  ind=obj\$ind,                                                #! was ist das?

call    = match.call(),
cov.mb  = value0\$cov.mb,
#                    cov.j   = cov.j,              # jackknife covariance matrix
#                    cov.lw  = cov.lw,             #
#                    cov.ls  = cov.ls,             #
#                    cov.method=cov.method,
#                    w.matrix= W\$w.matrix,         # weight matrix
#                    Wald    = value0\$outpar[9],
Wald =  (t(value0\$coefs) %*% solve(value0\$cov.rob)) %*% value0\$coefs,
means   = Means,               # means of <X> (model matrix)
linear.predictors= as.vector(scale(value0\$mmm, Means, scale=FALSE) %*% value0\$coefs),
method  = "Weighted Estimation",
method.ci= "Wald",             #
ci.lower= exp(value0\$coefs + qnorm(alpha/2) * vars^0.5), #
ci.upper= exp(value0\$coefs + qnorm(1 - alpha/2) * vars^0.5), #
prob    = probs,               # p-values
G       = G,
#                    W       = W,
W       = cbind(obj\$stratum,obj\$resp[,2],W)[obj\$resp[,3]==1,],
offset.values= obj\$offset.values, #
x       = if(x) obj\$mm1 else NA,   # return original model matrix if requested
npar = npar
)

#        ## if all weights are the same ...
#        if(W\$const) {
#                fit\$loglik <- value0\$outpar[12:11]
#                fit\$score <- value0\$outpar[7]
#        }

##        if(offset) {
##          obj\$covnames[1]<-as.character(paste("offset(",obj\$covnames[1],")"))
##          fit\$Wald <- t(coefs[2:k2]) %*% solve(covs[2:(k2),2:(k2)]) %*% coefs[2:(k2)]
##          if(x) {
##           fit\$x <- mm1.orig[,2:k2]
##           fit\$offset <- mm1.orig[,1]
##          }
##        }
names(fit\$prob) <- names(fit\$ci.upper) <- names(fit\$ci.lower) <- obj\$covnames
attr(fit, "class") <- c("concreg")
fit
}

coef.concreg<-function(object, ...) object\$coefficients

vcov.concreg<-function(object, ...) object\$var

cindex<-function(object, confint=FALSE, level=0.95) {
OC<-exp(coef(object))
cindex<-OC/(1+OC)
if(confint) {
ci<-confint(object, level=level, what="cindex")
cindex<-rbind(cindex,ci)
}
return(cindex)
}

confint.concreg<-function(object, parm, level=0.95, what="coefficients", ...){
if(what != "coefficients" & what !="OC" & what !="cindex") {
what<-"coefficients"
warning("Will produce confidence intervals for regression coefficients.\n")
}
if(missing(parm)) parm<-1:length(coef(object))
se<-sqrt(diag(vcov(object)))
ci.coef<-(coef(object)+se*cbind(qnorm((1-level)/2), qnorm(1-(1-level)/2))  )[parm,]
if(what=="coefficients") return(ci.coef)
else if (what=="OC") return (exp(ci.coef))
else if (what=="cindex") return (exp(ci.coef)/(1+exp(ci.coef)))
}

"plotw" <- function                                                             #! noch anpassen
(
x,    # object of class coxphw
rank=FALSE,
log=FALSE,
xlim=c(0,max(time)),

...            # dummy
)
{
if(rank) {
time<-order(x\$W[,2])
label<-"Ranked time"
}
else {
time<-x\$W[,2]
label<-"Time"
}
if(log) {
weights<-log(x\$W[,3])
wlabel<-"Log of weight"
}
else {
weights<-x\$W[,3]
wlabel<-"Weight"
}
#   if (is.na(xlim)) xlim=c(0,max(time))
ltyi<-1
for (i in unique(x\$W[,1])) {
if (ltyi==1) plot(time[x\$W[,1]==i], weights[x\$W[,1]==i],type="l",lty=ltyi,ylim=c(min(weights[time<=xlim[2]]),max(weights[time<=xlim[2]])), xlab=label, ylab=wlabel, xlim=xlim)
else lines(time[x\$W[,1]==i], weights[x\$W[,1]==i],lty=ltyi,ylim=c(min(weights[time<=xlim[2]]),max(weights[time<=xlim[2]])), xlab=label, ylab=wlabel, xlim=xlim)
ltyi<-ltyi+1
}
#   lines(time, weights,lty=2)
#   lines(time, weights,lty=3)
legend(min(time),0.95*max(weights[time<=xlim[2]]),unique(x\$W[,1]),lty=1:(ltyi-1), lwd=1)
}

"print.concreg" <- function                                                   #! noch anpassen
(
x,     # object of class coxphw
...            # dummy
)
### MP and GH
### overall test is score test
### 2007-07
{
print(x\$call)
cat("Model fitted by", x\$method, "\n\n")
if(x\$npar) cat("Nonparametric estimation \n")
se<- diag(x\$var)^0.5
out <- cbind(x\$coefficients, se, exp(x\$coefficients),
x\$ci.lower, x\$ci.upper, x\$coefficients/se, x\$prob)
dimnames(out) <- list(names(x\$coefficients),
c("coef", "se(coef)", "exp(coef)",
paste(c("lower", "upper"), 1 - x\$alpha), "z", "p"))

if (x\$method.ci != "Wald")
dimnames(out)[[2]][6] <- "Chisq"
print(out)

#        if("loglik" %in% names(x)) cat("\nScore test=", x\$score, " on ", x\$df,
#            " df, p=", 1 - pchisq(x\$score, x\$df), ", n=", x\$n, "", sep = "")
cat("\nWald test=", x\$Wald, " on ", x\$df, "df, p=", 1 - pchisq(x\$Wald, x\$df), ", n=", x\$n, "\n\n", sep = "")

invisible(x)
}

"summary.concreg" <- function                                                 #! noch anpassen
(
object,              # object of class coxphf
...                  # dummy
)
### MP and GH
### 2007-07
{
print(object\$call)
cat("\nModel fitted by", object\$method, "\n\n")
if(object\$npar) cat("Nonparametric estimation \n")
##cat("Confidence intervals and p-values by", object\$method.ci, "\n\n")
se<-diag(object\$var)^0.5
c.index<-exp(object\$coefficients)/(1+exp(object\$coefficients))
c.index.lower<-object\$ci.lower/(1+object\$ci.lower)
c.index.upper<-object\$ci.upper/(1+object\$ci.upper)
out <- cbind(object\$coefficients, se,
exp(object\$coefficients),
object\$ci.lower, object\$ci.upper,
c.index, c.index.lower, c.index.upper, object\$coefficients/se, object\$prob)
dimnames(out) <- list(names(object\$coefficients),
c("coef", "se(coef)", "exp(coef)",
paste(c("lower", "upper"), 1 - object\$alpha), "c", "c (lower)", "c (upper)", "z", "p"))
if (object\$method.ci != "Wald")
dimnames(out)[[2]][6] <- "Chisq"

print(out)

if("loglik" %in% names(object)) {
LL <- 2 * diff(object\$loglik)
cat("Likelihood ratio test=", LL, " on ", object\$df,
" df, p=", 1 - pchisq(LL, object\$df), ", n=", object\$n, "\n", sep = "")
cat("\nScore test=", object\$score, " on ", object\$df,
" df, p=", 1 - pchisq(object\$score, object\$df), ", n=", object\$n, "\n", sep = "")
}
#        wald.z <- t(coef(object)) %*% solve(object\$var) %*% coef(object)
cat("Wald test =", object\$Wald, "on", object\$df, " df, p =",
1 - pchisq(object\$Wald, object\$df))
cat("\n\nCovariance-Matrix:\n")
print(object\$var)

invisible(object)
}

concreg.fit <- function(obj, id, W, G, PARMS, npar)
{
# fitter function

#        k <- ncol(obj\$mm1)    # number covariates w/o time-dep effects         #! in PARMS
k <- PARMS[2]
k2 <- k + obj\$NTDE                                                      #! weg?
#        maxid <- max(id)                                                       #! in PARMS
maxid <- PARMS[15]

## standardize model matrix, but only in semiparametric mode
if(!npar) {
sd1 <- apply(as.matrix(obj\$mm1),2,sd)
sd2 <- apply(as.matrix(obj\$timedata),2,sd)                                                 #! weg?
Z.sd <- c(sd1, sd2 * sd1[obj\$timeind])                                  #! weg?
ZxZ <- as.matrix(Z.sd) %*% t(as.matrix(Z.sd)) # used often to restandardize ...    #! weg?
obj\$mm1 <- scale(obj\$mm1, FALSE, sd1)
} else {
sd1<-1
sd2 <- 1
Z.sd<-1
ZxZ <-1
}

##if(ind.offset)
obj\$mm1o <- if(PARMS[16] != 0) cbind(obj\$offset.values, obj\$mm1) else obj\$mm1

#        if(is.null(CARDS))
CARDS <- cbind(obj\$mm1o, obj\$resp[,c(2, 4)], W, G, id, obj\$stratum)
#          CARDS <- cbind(obj\$mm1o, obj\$resp, weights, obj\$timedata, id)        #! was ist weights?

if(!npar) obj\$timedata <- scale(obj\$timedata, FALSE, sd2)                         # weg?
mmm <- cbind(obj\$mm1, obj\$timedata) # model matrix inc. time data       #! ohne timedata?
##   if (offset) {
##    IOARRAY[1,1]<-0    # first variable is offset
##    IOARRAY[2,1]<-Z.sd[1]    # first variable is offset
##   }
DFBETA <- matrix(0, maxid, k2)                                          #! k2 durch k ersetzen?
IOARRAY <- rbind(rep(1, k2), matrix(0, 2+2*k2, k2))                     #! k2 durch k ersetzen?
if(obj\$NTDE >0)                                                         #! weg?
IOARRAY[4, (k+1):k2] <- obj\$timeind                                   #! weg?

## --------------- Aufruf Fortran-Routine ????????? ------------
storage.mode(CARDS) <- storage.mode(PARMS) <- storage.mode(IOARRAY) <- storage.mode(DFBETA) <- "double"
value <- .Fortran("CONCREG",                                            #! anpassen
cards=CARDS,
outpar = PARMS,
outtab = IOARRAY)
#                          PACKAGE=concreg.fp)
if(value\$outpar[8])
warning("Error in Fortran routine concreg; parms8 <> 0")

coefs <- value\$coefs / Z.sd                                             #! value\$coefs gibt es nicht
cov.mb <- matrix(value\$outtab[4:(k2+3), ], ncol=k2) / ZxZ            # "model-based" varianz
cov.rob <- matrix(value\$outtab[(3+k2+1):(3+2*k2), ], ncol=k2) / ZxZ  # "robuste" varianz (standard)

res <- list(                                                            #! anpassen
cards=value\$cards,
outpar=value\$outpar,
outtab=matrix(value\$outtab, nrow=3+2*k2),
#                   dfbetaresid=value\$dfbetaresid,
coef.orig=value\$outtab[3,  ],
coefs=value\$outtab[3,  ] / Z.sd, # coefficients
cov.rob=cov.rob,                 # covariances
cov.mb=cov.mb,
Z.sd=Z.sd,
ZxZ=ZxZ,                                                    #! weg?
mmm=mmm                          # model matrix
)
res
}

concreg.wei <- function(resp, max.strata, stratum, trunc.weights)
{
#! function to calculate the weigth for concreg
#! rechnet mit time.crisk

data.tmp <- data.frame(resp[,c(5,3,2,4)])
dimnames(data.tmp)[[2]][2] <- "status"
dimnames(data.tmp)[[2]][3] <- "time.orig"
dimnames(data.tmp)[[2]][4] <- "status012"
data.tmp\$statusfu<-(data.tmp\$status==0)+0

W <- rep(0, nrow(data.tmp))
if(all(data.tmp[,2]==1)) return(rep(1, nrow(data.tmp)))

for (i.max.strata in 1:max.strata)  {
data.strata <- data.tmp[stratum==i.max.strata,]

n.strata <- nrow(data.strata)
data.strata\$pat <- 1:n.strata

#    attach(data.strata, warn.conflicts=FALSE)         # survfit geht sonst nicht aus irgendeinem Grund
fit <- survfit(Surv(time.crisk, status)~1, data=data.strata)           #pseudo survival to reconstruct number of pairs
fit.fu<-survfit(Surv(time.orig, statusfu)~1, data=data.strata)         #follow-up KM
#    detach(data.strata)
fsurv<-fit\$surv

#   if (length(fsurv) == n.strata)  {                        # keine Bindungen
#      fsurv <- c(1, fit\$surv)[-(length(fit\$surv)+1)]
#     data.strata\$surv <- fsurv
#      data.strata\$n.risk <- fit\$n.risk
#    }

#   if (length(fsurv) != n.strata)  {                        # Bei Bindungen innerhalb der Strata
data.strata\$surv <- -99
data.strata\$n.risk <- -99
data.strata\$gsurv <- -99

fit.n <- 1
fsurv<-c(1,fit\$surv)
fn.risk<-c(n.strata, fit\$n.risk)
gsurv<-c(1,fit.fu\$surv)
gtime<-c(0,fit.fu\$time)
ftime<-c(0,fit\$time)
for (i.n in 1:n.strata)  {
if (data.strata\$status[i.n]==1)  indices<-sum(ftime<data.strata[i.n, "time.crisk"])
if (data.strata\$status[i.n]==0)  indices<-sum(ftime<=data.strata[i.n, "time.crisk"])
index.fu<-sum(gtime<data.strata[i.n, "time.crisk"])
data.strata[i.n, c("surv", "n.risk", "gsurv")] <- c((fsurv[indices]), (fn.risk[indices]), gsurv[index.fu])
#        if (data.strata[i.n, "time.crisk"] == fit\$time[fit.n]) {
#          fit.n <- fit.n+1
}
#      }

#    data.strata[data.strata[order(data.strata[,"surv"])[1:sum(data.strata[,"surv"] == -99)],3], c("surv", "n.risk")] <-
#    data.strata[data.strata[order(data.strata[,"surv"])[1:sum(data.strata[,"surv"] == -99)]-1,3], c("surv", "n.risk")]

#    W <- c(W, (n.strata*data.strata[,"surv"] / (data.strata[,"n.risk"]-1))* (data.strata[,"surv"] / (data.strata[,"n.risk"])))
#    W <- c(W, ((n.strata*data.strata[,"surv"]-1) / (data.strata[,"n.risk"]-1))* (n.strata*data.strata[,"surv"] / (data.strata[,"n.risk"])))
W.strata <- ((n.strata*data.strata[,"surv"]-1) / (data.strata[,"n.risk"]-1))/ data.strata[,"gsurv"]
W[stratum==i.max.strata] <- W.strata

}
#  W <- W[-1]
W[W==Inf] <- 1
W[is.na(W)] <- 1
truncat<-quantile(W, probs=trunc.weights)
W[W>truncat]<-truncat
W[abs(W)==Inf]<-0
W
}
```

## Try the concreg package in your browser

Any scripts or data that you put into this service are public.

concreg documentation built on May 29, 2017, 1:41 p.m.