0

我需要创建自己的 nribin 代码，它可以用于 logistf 包结果，我也许可以工作，请给我一些建议！！！

我将 z.std = mdl.std$x[,-1] 更改为 z.std = mdl.std$x 并取消：link = mdl.std$family[[2]] 和 family=binomial(link)，

整个代码是：

nribin_LTY <- 函数（事件=NULL，mdl.std=NULL，mdl.new=NULL，z.std=NULL，z.new=NULL，p.std=NULL，p.new=NULL，updown='category' ，cut=NULL，link='logit'，niter=1000，alpha=0.05，msg=TRUE）{

##
## type of calculation
flag.mdl = !is.null(mdl.std) && !is.null(mdl.new)
flag.prd = !is.null(z.std)   && !is.null(z.new)
flag.rsk = !is.null(p.std)   && !is.null(p.new)



##
## check standard & new model
if (flag.mdl) {
  if (is.null(event)) event = as.numeric(mdl.std$y)
  if (is.null(mdl.std$x) || is.null(mdl.new$x))
    stop("\n\nmodel object does not contain predictors. pls set x=TRUE for model calculation.\n\n")

  z.std   = mdl.std$x
  z.new   = mdl.new$x

  mdl.std = glm(event ~ .,  data=as.data.frame(cbind(event, z.std)))
  mdl.new = glm(event ~ .,  data=as.data.frame(cbind(event, z.new)))

} else if (flag.prd) {
  mdl.std = glm(event ~ .,  data=as.data.frame(cbind(event, z.std)))
  mdl.new = glm(event ~ .,  data=as.data.frame(cbind(event, z.new)))***
  message("\nSTANDARD prediction model:")
  print(summary(mdl.std)$coef)
  message("\nNEW prediction model:")
  print(summary(mdl.new)$coef)

} else if (!flag.mdl && !flag.prd && !flag.rsk) {
  stop("\n\neither one of 'event, z.std, z.new', 'event, p.std, p.new', and 'mdl.std, mdl.new' should be specified.\n\n")
}

if (is.null(cut))
  stop("\n\n'cut' is empty")

objs = list(mdl.std, mdl.new, z.std, z.new, p.std, p.new)

##
## DH & DL
wk  = get.uppdwn.bin(event, objs, flag.mdl, flag.prd, flag.rsk, updown, cut, link, msg=msg)
upp = wk[[1]]
dwn = wk[[2]]
ret = list(mdl.std=mdl.std, mdl.new=mdl.new, p.std=wk[[3]], p.new=wk[[4]], up=upp, down=dwn, rtab=wk[[5]], rtab.case=wk[[6]], rtab.ctrl=wk[[7]])

##
## point estimation
message("\nNRI estimation:")
est = nribin.count.main(event, upp, dwn)
message("Point estimates:")
result = data.frame(est)
names(result) = 'Estimate'
row.names(result) = c('NRI','NRI+','NRI-','Pr(Up|Case)','Pr(Down|Case)','Pr(Down|Ctrl)','Pr(Up|Ctrl)')
print(result)

##
## interval estimation
if (niter > 0) {
  message("\nNow in bootstrap..")
  ci   = rep(NA, 14)
  N    = length(event)
  samp = matrix(NA, niter, 7)
  colnames(samp) = c('NRI','NRI+','NRI-','Pr(Up|Case)','Pr(Down|Case)','Pr(Down|Ctrl)','Pr(Up|Ctrl)')

  for (b in 1:niter) {
    f    = as.integer(runif(N, 0, N)) + 1
    objs = list(mdl.std, mdl.new, z.std[f,], z.new[f,], p.std[f], p.new[f])
    wk   = get.uppdwn.bin(event[f], objs, flag.mdl, flag.prd, flag.rsk, updown, cut, link, msg=FALSE)
    upp  = wk[[1]]
    dwn  = wk[[2]]
    samp[b,] = nribin.count.main(event[f], upp, dwn)
  }
  ret = c(ret, list(bootstrapsample=samp))

  ci = as.numeric(apply(samp, 2, quantile, c(alpha/2, 1-alpha/2), na.rm=TRUE, type=2))
  se = as.numeric(apply(samp, 2, sd))
  message("\nPoint & Interval estimates:")
  result = as.data.frame(cbind(est, se, matrix(ci, ncol=2, byrow=TRUE)))
  names(result) = c('Estimate', 'Std.Error', 'Lower', 'Upper')
  row.names(result) = c('NRI','NRI+','NRI-','Pr(Up|Case)','Pr(Down|Case)','Pr(Down|Ctrl)','Pr(Up|Ctrl)')
  print(result)
}

invisible(c(list(nri=result), ret))

}

4

1 回答 1

0

需要在logistf结果中添加变量矩阵，例如：

x.std=as.matrix(thromb2[,c(8,9,10,14,21,22,24,25,26)])

mstd = logistf(formula=fml.std, firth = FALSE,
               data = thromb2, x=TRUE)
mstd =backward(mstd)
mstd$x <- x.std
于 2018-01-21T14:28:02.797 回答