Nothing
# type 3 and 7 do not give the right output for glm fits
# robust can be passed as part of .... Sometimes robust=T generates errors
# trim: get rid of white space in confidence intervals for alignment
getFormattedSummary=function(fits, type=12, est.digits=2, se.digits=2, robust, random=FALSE, VE=FALSE, to.trim=FALSE, rows=NULL, coef.direct=FALSE, trunc.large.est=TRUE, scale.factor=1, p.digits=3, remove.leading0=FALSE, p.adj.method="fdr", ...){
if(is.null(names(fits))) names(fits)=seq_along(fits)
idxes=seq_along(fits); names(idxes)=names(fits)
if(length(robust)==1) robust=rep(robust,length(fits)) else if (length(robust)!=length(fits)) stop("length of robust needs to match length of fits")
res = sapply(idxes, simplify="array", function (fit.idx) {
fit=fits[[fit.idx]]
if(coef.direct) {
tmp=fit
} else {
if (random) {
tmp = getVarComponent (fit)
if (inherits(fit, "mer") & type!=1) {
warning ("only point estimate is available for variance components from lmer fit, forcing type to 1")
type=1
}
} else {
tmp = getFixedEf (fit, robust=robust[fit.idx], ...)
}
if (VE) {
tmp[,1]=1-tmp[,1]
# reverse lb and ub
tmpv=tmp[,4]
tmp[,4]=1-tmp[,3]
tmp[,3]=1-tmpv
}
}
# tmp should be: est, se, lb, ub, pvalue
tmp[,1]=tmp[,1]*scale.factor
if(ncol(tmp)>1) tmp[,2]=tmp[,2]*scale.factor
# take only some rows
if (!is.null(rows)) tmp=tmp[intersect(rows,1:nrow(tmp)),,drop=F]
est.=as.matrix(tmp[,1,drop=FALSE], ncol=1) # if do not cast into matrix, it will be a data frame, which leads to problems next
too.big=trunc.large.est & est.>100
# replace large values in est.
# find a round that takes multipled digits!
# trim is necessary here
est.=ifelse(too.big,">100",trim(formatDouble(est., est.digits, remove.leading0=remove.leading0)))
p.val.col=which(startsWith(tolower(colnames(tmp)),"p"))
if(ncol(tmp)>1) {
lb=tmp[,3,drop=FALSE]
lb[too.big]=NA
lb=formatDouble(lb, est.digits, remove.leading0=remove.leading0)
if(to.trim) lb=trim(lb)
ub=tmp[,4,drop=FALSE]
ub[too.big]=NA
ub=formatDouble(ub, est.digits, remove.leading0=remove.leading0)
if(to.trim) ub=trim(ub)
}
# str(lb); str(ub); str(est.) # make sure they are not data frames
if (type==1)
# est only
out=drop(est. )
else if (type==2)
# est (se)
out=est. %.% " (" %.% formatDouble(tmp[,2,drop=FALSE], se.digits, remove.leading0=remove.leading0) %.% ")" %.% ifelse (round(tmp[,p.val.col],3)<=0.05, ifelse (tmp[,p.val.col]<0.01,"**","*"),"")
else if (type==3)
# est (lb, up)
out=est. %.% " (" %.% lb %.% ", " %.% ub %.% ")"
else if (type==4)
# a space is inserted between est and se, they could be post-processed in swp
out=est. %.% " " %.% formatDouble(tmp[,2,drop=FALSE], est.digits, remove.leading0=remove.leading0)
else if (type==5)
# est **
out=est. %.%
ifelse (round(tmp[,p.val.col],3)<=0.05,ifelse (tmp[,p.val.col]<0.01,"**","*"),"")
else if (type==6)
# est (pval)*
out=est. %.% " (" %.% formatDouble(tmp[,p.val.col,drop=FALSE], 3, remove.leading0=remove.leading0) %.% ")" %.% ifelse (round(tmp[,p.val.col],3)<=0.05,ifelse (tmp[,p.val.col]<0.01,"**","*"),"")
else if (type==7)
# (lb, up)
out=ifelse(drop(too.big), rep("",nrow(lb)), " (" %.% lb %.% ", " %.% ub %.% ")")
else if (type==8)
# est (p value #)
out=est. %.% " (p value " %.%
formatDouble(tmp[,p.val.col,drop=FALSE], 3, remove.leading0=remove.leading0) %.% ")"
else if (type==9)
# est (pval)*
out=est. %.% " (" %.% formatDouble(tmp[,p.val.col,drop=FALSE], p.digits, remove.leading0=remove.leading0) %.% ")"
else if (type==10) {
# pval
tmp.out=tmp[,p.val.col,drop=TRUE]
out = ifelse(tmp.out<10^(-p.digits), paste0("<0.",concatList(rep("0",p.digits-1)),"1"), format(round(tmp.out, p.digits), nsmall=p.digits, scientific=FALSE) )
}
else if (type==11)
# adj pval
out=format(round(p.adjust(tmp[,p.val.col,drop=TRUE], method=p.adj.method), p.digits), nsmall=3, scientific=FALSE)
else if (type==12)
# est (lb, up, pval *)
out=est. %.%
" (CI=" %.% lb %.% "," %.% ub %.%
", p=" %.% formatDouble(tmp[,p.val.col,drop=FALSE], 3, remove.leading0=remove.leading0) %.% ")" %.%
ifelse (round(tmp[,p.val.col],3)<=0.05,ifelse (tmp[,p.val.col]<0.01,"**","*"),"")
else if (type==13)
# (lb-up)
out="(" %.% lb %.% "-" %.% ub %.% ")"
else
stop ("getFormattedSummaries(). type not supported: "%.%type)
names(out)=rownames(tmp)
out=gsub("NA","",out)
out=gsub("\\( +\\)","",out)
out
})
if (is.list(res)) {
# if the fits have different number of parameters, we need this
res=cbinduneven(li=lapply(res, function (x) as.matrix(x, ncol=1)))
colnames(res)=names(fits)
} else if (!is.matrix(res)) {
# if there is only one coefficient, we need this
res=matrix(res, nrow=1)
colnames(res)=names(fits)
}
res
}
# return a matrix, columns are: est, p value, low bound, upper bound
getFixedEf <- function(object, ...) UseMethod("getFixedEf")
# variance component
getVarComponent <- function(object, ...) UseMethod("getVarComponent")
# if there is missing data and robust=T, some errors will happen. it is a useful error to have.
# if ret.robcov TRUE, then returns robust variance-covariance matrix
# infjack.glm defined later in this file
getFixedEf.glm = function (object, exp=FALSE, robust=TRUE, ret.robcov=FALSE, ...) {
x=summary(object)
out=x$coef
if (robust | ret.robcov) {
V=infjack.glm(object, 1:nrow(object$model)) # object$data may have NAs
if (ret.robcov) return (V)
out[,2]=sqrt(diag(V))
out[,3]=out[,1]/out[,2]
out[,4]=2 * pnorm(-abs(out[,3]))
}
# to get est, p value, low bound, upper bound
out=cbind(out[,1],out[,2],out[,1]-1.96*out[,2],out[,1]+1.96*out[,2],out[,4])
colnames(out)=c("est","se","(lower","upper)","p.value")
if(exp) {
out[,c(1,3,4)]=exp(out[,c(1,3,4)])
colnames(out)[1]="OR"
}
# sometimes some coef are null, we may want to preserve, modified from print.summary.glm
if(!is.null(aliased <- x$aliased) && any(aliased)) {
cn <- names(aliased)
coefs <- matrix(NA, length(aliased), ncol(out),
dimnames=list(cn, colnames(out)))
coefs[!aliased, ] <- out
out = coefs
}
out
}
getFixedEf.gee = function (object,exp=FALSE, ...) {
#print("in getFixedEf.gee")
out=as.matrix(summary(object)$coef )
# Estimate Std.err Wald Pr(>|W|)
out=cbind(out[,1:2], out[,1]-1.96*out[,2], out[,1]+1.96*out[,2], out[,4,drop=FALSE])
colnames(out)=c("est","se","(lower","upper)","p.value")
if(exp) {
out[,c(1,3,4)]=exp(out[,c(1,3,4)])
colnames(out)[1]="RR"
}
out
}
getFixedEf.MIresult=function(object, ...) {
capture.output({
tmp=summary(object)
}, type="output") # type = message captures stderr, type=output is for stdout
tmp=tmp[,names(tmp)!="missInfo"] # MIresult has this string column #tmp=subset(tmp, select=-missInfo) fails check
tmp=as.matrix(tmp)# data frame fails in format and formatDouble
cbind(tmp, "p.val"=2*pt(abs(tmp[,1]/tmp[,2]), df=object$df, lower.tail = FALSE))
}
#get estimates, variances, sd from lmer fit
getFixedEf.mer = function (object, ...) {
Vcov <- lme4::VarCorr(object, useScale = FALSE)
betas <- lme4::fixef(object)
se <- sqrt(diag(Vcov))
zval <- betas / se
pval <- 2 * pnorm(abs(zval), lower.tail = FALSE)
cbind("Estimate"=betas, se, zval, pval)
}
getVarComponent.mer = function (object, ...) {
tmp=lme4::VarCorr(object)
mysapply(tmp, function (comp) attr(comp, "stddev") )
}
#get estimates, variances, sd from lmer fit
getFixedEf.glmerMod = function (object, exp=FALSE, ...) {
out=summary(object)$coefficients
out=cbind(est=out[,1], se=out[,2], lb=out[,1]-1.96*out[,2], ub=out[,1]+1.96*out[,2],"p"=out[,4])
if(exp) out[,c(1,3,4)]=exp(out[,c(1,3,4)])
out #est, se, lb, ub, pvalue
}
getFixedEf.merMod=getFixedEf.glmerMod
getVarComponent.glmerMod = function (object, ...) {
tmp=lme4::VarCorr(object)
mysapply(tmp, function (comp) attr(comp, "stddev") )
}
getVarComponent.merMod=getVarComponent.glmerMod
#get estimates, variances, sd from lme fit
getFixedEf.lme = function (object, ...) {
betas <- object$coef$fixed
se <- sqrt (diag (object$varFix))
zval <- betas / se
pval <- 2 * pnorm(abs(zval), lower.tail = FALSE)
cbind(betas, se, zval, pval)
}
# getVarComponent.lme = function (object, ...) {
# nlme::VarCorr(object)
# }
getFixedEf.lmerMod = function (object, ...) {
betas <- nlme::fixef(object)
se <- sqrt (diag (getVarComponent(object)))
zval <- betas / se
pval <- 2 * pnorm(abs(zval), lower.tail = FALSE)
cbind(betas, se, zval, pval)
}
getVarComponent.lmerMod = function (object, ...) {
as.matrix(vcov(object)) # otherwise will complain about S4 class convertibility problem
}
getFixedEf.geese = function (object, robust=TRUE, ...) {
tmp=summary(object)$mean
# tmp is: estimate san.se wald p
# tmp should be: est, se, lb, ub, pvalue
as.matrix(cbind(tmp[,c(1, ifelse(robust,2,3))], confint(object), tmp[,4,drop=F])) # if leave as data.frame, formatDouble fails
}
getFixedEf.logistf = function (object, exp=FALSE, ...) {
temp = summary(object)
out = cbind (coef=temp$coef, se=NA, p.value=temp$prob)
if(exp) out[,1]=exp(out[,1])
out
}
vcov.logistf=function(object, ...){
object$var
}
getFixedEf.gam = function (object, ...) {
temp = summary(object)
cbind (temp$p.coef, temp$se[1:length(temp$p.coef)])
}
getFixedEf.lm = function (object, ...) {
out=summary(object)$coef
ci=confint(object)
out=cbind(out[,1:2], ci, out[,4])
colnames(out)[5]="p-val"
out
}
getFixedEf.inla = function (object, ...) {
tmp = summary(object)$fixed
n=nrow(tmp)
tmp.name = row.names(tmp)[n]
# move intercept to the top
if (tmp.name=="intercept") {
tmp = rbind (tmp[n,],tmp)[1:n,,drop=FALSE]
dimnames (tmp)[[1]][1] = tmp.name
}
# rename first column
dimnames (tmp)[[2]][1] = "Estimate"
tmp
}
# return the mean, sd, CI of the transformed variable
inla.getMeanSd=function (marginal, f="identity") {
# interpolations suggested by Havard: do it on the original scale
logtao=log(marginal[,1]); p.logtao=marginal[,2]*marginal[,1]
fun = splinefun(logtao, log(p.logtao))
h=0.001
x = seq(min(logtao),max(logtao),by=h)
pmf = exp(fun(x))*h
# sum (pmf) # Prob
# x = seq(min(logtao)-sd(logtao)/2,max(logtao)+sd(logtao)/2,by=h)
# pmf = exp(fun(x))*h
# sum (pmf) # Prob
# x = seq(min(logtao)-sd(logtao),max(logtao)+sd(logtao),by=h)
# pmf = exp(fun(x))*h
# sum (pmf) # Prob
# x = seq(min(logtao)-sd(logtao)*2,max(logtao)+sd(logtao)*2,by=h)
# pmf = exp(fun(x))*h
# sum (pmf) # Prob
lower.boundary = rle(cumsum(pmf)>.025)$lengths[1]+1
upper.boundary = rle(cumsum(pmf)>.975)$lengths[1]+1
# if (pmf[lower.boundary]>.04) {
# #stop ("getMeanSd(): pmf too large at lower boundary: "%.%pmf[lower.boundary])
# return (rep(NA, 4))
# }
# if (pmf[upper.boundary]>.04) {
# stop ("getMeanSd(): pmf too large at upper boundary"%.%pmf[upper.boundary])
# return (rep(NA, 4))
# }
if (f=="identity") {
func=function(x) { exp(x) }
} else if (f=="inverse") {
func=function(x) { exp(x)**-1 }
} else if (f=="inversesqrt") {
func=function(x) { exp(x)**-.5 }
} else
stop ("getMeanSd(): function not supported "%.%f)
mu = sum( pmf * func(x))
stdev = (sum( pmf * func(x)**2 ) - mu**2) ** .5
out = c("mean"=mu, "stddev"=stdev, 0,0 ) # we may run into problems with the boundary
#out = c("mean"=mu, "stddev"=stdev, sort(func(x)[c(lower.boundary, upper.boundary)]) );
names(out)[3]="2.5%"; names(out)[4]="97.5%";
out
}
# returns estimate of standard deviation and the estimated sd of that estimate
getVarComponent.hyperpar.inla = function (object, transformation=NULL, ...) {
marginals = object$marginals
out = mysapply(1:length(marginals), function (i) {
# this is a little precarious, but hey
if (startsWith(names(marginals)[i],"Prec")) {
if (is.null (transformation)) {
inla.getMeanSd(marginals[[i]],"inversesqrt")
} else {
inla.getMeanSd(marginals[[i]],transformation)
}
} else if (startsWith(names(marginals)[i],"Rho")) {
object$summary[i, c(1,2,3,5)]
} else {
stop ("don't know what to do with this names(marginals)[i]: "%.% names(marginals)[i] )
}
})
dimnames (out)[[1]]="sigma.or.rho."%.%dimnames (out)[[1]]
out
}
getFixedEf.svycoxph=function (object, exp=FALSE, robust=TRUE, ...){
if (!robust) warning("svycoxph variance estimate is always design-based, which is close to robust variance estimate. No model-based variance estimate is available")
hr=object$coefficients
se=sqrt(diag(object$var))
pval=pnorm(abs(hr/se), lower.tail=F)*2
out=cbind(HR=hr,
"se"=se,
"(lower"=hr-qnorm(0.975)*se,
"upper)"=hr+qnorm(0.975)*se,
"p.value"=pval
)
if (exp) out[,c(1,3,4)]=exp(out[,c(1,3,4)])
out
# round(sqrt(diag(attr(object$var,"phases")$phase1)),3)
# round(sqrt(diag(attr(object$var,"phases")$phase2)),3)
}
getFixedEf.svyglm=function (object, exp=FALSE, robust=TRUE, ...){
if (!robust) warning("svyglm variance estimate is always design-based, which is close to robust variance estimate. No model-based variance estimate is available")
tmp=summary(object)$coefficients
or=tmp[,"Estimate"]
se=tmp[,"Std. Error"]
pval=tmp[,"Pr(>|t|)"]
out=cbind(OR=or,
"se"=se,
"(lower"=or-qnorm(0.975)*se,
"upper)"=or+qnorm(0.975)*se,
"p.value"=pval
)
if (exp) out[,c(1,3,4)]=exp(out[,c(1,3,4)])
out
# round(sqrt(diag(attr(object$var,"phases")$phase1)),3)
# round(sqrt(diag(attr(object$var,"phases")$phase2)),3)
}
getFixedEf.svy_vglm=function (object, exp=FALSE, robust=TRUE, ...){
if (!robust) warning("svy_vglm variance estimate is always design-based, which is close to robust variance estimate. No model-based variance estimate is available")
tmp=summary(object)$coeftable
or=tmp[,"Coef"]
se=tmp[,"SE"]
pval=tmp[,"p"]
out=cbind(HR=or,
"se"=se,
"(lower"=or-qnorm(0.975)*se,
"upper)"=or+qnorm(0.975)*se,
"p.value"=pval
)
if (exp) out[,c(1,3,4)]=exp(out[,c(1,3,4)])
out
}
getFixedEf.coxph=function (object, exp=FALSE, robust=FALSE, ...){
#capture.output()# summary.svycoxph prints some stuff, capture.output absorbs it
sum.fit<-summary(object)
se.idx=ifelse(!robust,"se(coef)","robust se")
if(robust) {
pvals=sum.fit$coef[,"Pr(>|z|)"]
} else {
pvals=pnorm(abs(sum.fit$coef[,1]/sum.fit$coef[,"se(coef)"]), lower.tail=F)*2
}
if (!exp) {
cbind(HR=sum.fit$coef[,1],
"se"=sum.fit$coef[,se.idx],
"(lower"=(sum.fit$coef[,1]-qnorm(0.975)*sum.fit$coef[,se.idx]),
"upper)"=(sum.fit$coef[,1]+qnorm(0.975)*sum.fit$coef[,se.idx]),
"p.value"=pvals
)
} else {
cbind(HR=sum.fit$coef[,2],
"se"=sum.fit$coef[,se.idx],
"(lower"=exp(sum.fit$coef[,1]-qnorm(0.975)*sum.fit$coef[,se.idx]),
"upper)"=exp(sum.fit$coef[,1]+qnorm(0.975)*sum.fit$coef[,se.idx]),
"p.value"=pvals
)
}
# round(sqrt(diag(attr(object$var,"phases")$phase1)),3)
# round(sqrt(diag(attr(object$var,"phases")$phase2)),3)
}
getFixedEf.matrix = function (object, ...) {
t(apply(object, 2, function (x) c("Estimate"=mean(x), "sd"=sd(x), "2.5%"=quantile(x,.025), "97.5"=quantile(x,.975))))
}
getVarComponent.matrix = function (object, ...) {
t(apply(object, 2, function (x) c("Estimate"=mean(x), "sd"=sd(x), "2.5%"=quantile(x,.025), "97.5"=quantile(x,.975))))
}
#################################################################################################
coef.geese = function (object, ...) {
tmp = summary(object)$mean[,1]
names (tmp)=names (object$beta)
tmp
}
vcov.geese = function (object, ...) {
tmp = object$vbeta
dimnames (tmp)=list (names (object$beta), names (object$beta))
tmp
}
residuals.geese = function (object, y, x, ...) {
y - x %*% object$beta
}
predict.geese = function (object, x, ...) {
x %*% object$beta
}
#################################################################################################
coef.tps = function (object, ...) {
object$coef
}
vcov.tps = function (object, robust, ...) {
if(robust) object$cove else object$covm
}
getFixedEf.tps = function (object, exp=FALSE, robust=TRUE, ...) {
res = summary(object)$coef
idx=ifelse(robust,"Emp ","Mod ")
res = cbind(res[,c("Value",idx%.%"SE")], "lower bound"=res[,"Value"]-1.96*res[,idx%.%"SE"], "upper bound"=res[,1]+1.96*res[,idx%.%"SE"], "p value"=res[,idx%.%"p"])
if (exp) res[,c(1,3,4)] = exp(res[,c(1,3,4)])
colnames(res)=c(ifelse(exp,"OR","est"), "se(est)", "(lower", "upper)", "p.value")
res
}
predict.tps = function (object, newdata = NULL, type = c("link", "response"), ...) {
type=match.arg(type)
out=as.matrix(newdata) %*% coef.tps(object)
if(type=="response") out=expit(out)
out
}
# make two tables, rotating v1 and v2
# fit is an object that needs to have coef and vcov
# list the effect of one variable at three levels of another variable, if data is in fit, use quantiles, otherwise -1,0,1
# v1.type and v2.type: continuous, binary, etc
interaction.table=function(fit, v1, v2, v1.type="continuous", v2.type="continuous", logistic.regression=TRUE){
coef.=coef(fit)
cov.=vcov(fit)
var.names=names(coef.)
v1.ind = match(v1, var.names)
v2.ind = match(v2, var.names)
itxn.ind = match(v1%.%":"%.%v2, var.names)
if(is.na(itxn.ind)) itxn.ind = match(v2%.%":"%.%v1, var.names)
if (any(is.na(c(v1.ind, v2.ind, itxn.ind)))) {
stop("v1, v2, or interaction not found in var.names")
}
ret=list()
for (i in 1:2) {
if (i==1) {
ind.1=v1.ind; type.1=v1.type
ind.2=v2.ind; type.2=v2.type
} else {
ind.1=v2.ind; type.1=v2.type
ind.2=v1.ind; type.2=v1.type
}
# three levels of v.2
if(is.null(fit$data)) {
lev=c(-1,0,1)
names(lev)=lev
} else {
# get quantiles
lev=quantile(fit$data[[var.names[ind.2]]], c(.25,.5,.75), na.rm=TRUE)
}
# increase ind.1 by 1 at lev
lin.combs=NULL
if(type.2=="binary") {
lin.comb.1=rep(0, length(coef.)); lin.comb.1[ind.1]=1; lin.comb.1[itxn.ind]=0; lin.combs=rbind(lin.combs, "0"=lin.comb.1)
lin.comb.1=rep(0, length(coef.)); lin.comb.1[ind.1]=1; lin.comb.1[itxn.ind]=1; lin.combs=rbind(lin.combs, "1"=lin.comb.1)
} else {
lin.comb.1=rep(0, length(coef.)); lin.comb.1[ind.1]=1; lin.comb.1[itxn.ind]=lev[1]; lin.combs=rbind(lin.combs, lin.comb.1)
lin.comb.1=rep(0, length(coef.)); lin.comb.1[ind.1]=1; lin.comb.1[itxn.ind]=lev[2]; lin.combs=rbind(lin.combs, lin.comb.1)
lin.comb.1=rep(0, length(coef.)); lin.comb.1[ind.1]=1; lin.comb.1[itxn.ind]=lev[3]; lin.combs=rbind(lin.combs, lin.comb.1)
rownames(lin.combs)=names(lev)
}
effect = lin.combs%*%coef.
sd. = sqrt(diag(lin.combs%*%cov.%*%t(lin.combs)))
p.val = pnorm(abs(effect)/sd., lower.tail=FALSE)
lci=effect-1.96*sd.
uci=effect+1.96*sd.
res = cbind(effect, lci, uci, p.val)
colnames(res)=c("coef","(lower","upper)","p value")
if (logistic.regression) {
res[,1:3]=exp(res[,1:3])
colnames(res)[1]="OR"
}
ret[[i]]=res
}
names(ret) = "Effect of increasing "%.%c(v1,v2) %.% " by 1 at selected values of " %.% c(v2,v1)
ret
}
# these functions are needed by some of the functions in this file
# returns sandwich estimator of variance matrix
# from Thomas Lumley
infjack.glm<-function(glm.obj,groups){
umat<-estfun.glm(glm.obj)
usum<-rowsum(umat,groups,reorder=FALSE)
modelv<-summary(glm.obj)$cov.unscaled
modelv%*%(t(usum)%*%usum)%*%modelv
}
jack.glm<-function(glm.obj,groups){
umat<-jackvalues(glm.obj)
usum<-rowsum(umat,groups,reorder=FALSE)
t(usum)%*%usum*(nrow(umat)-1)/nrow(umat)
}
jackvalues<-function(glm.obj){
db<-lm.influence(glm.obj)$coef
t(t(db)-apply(db,2,mean))
}
estfun.glm<-function(glm.obj){
if (is.matrix(glm.obj$x))
xmat<-glm.obj$x
else {
mf<-model.frame(glm.obj)
xmat<-model.matrix(terms(glm.obj),mf)
}
residuals(glm.obj,"working")*glm.obj$weights*xmat[,rownames(summary(glm.obj)$cov.unscaled)] # the column selection is added to fix a bug caused by weird formula
}
# goodness of fit
# ngroups=10; main=""; add=FALSE; show.emp.risk=TRUE; lcol=1; ylim=NULL; weights=NULL
risk.cal=function(risk, binary.outcome, weights=NULL, ngroups=NULL, cuts=NULL, main="", add=FALSE, show.emp.risk=TRUE, lcol=2, ylim=NULL, scale=c("logit","risk")){
scale=match.arg(scale)
stopifnot(length(risk)==length(binary.outcome))
if(!is.null(weights)) stopifnot(length(risk)==length(weights)) else weights=rep(1,length(risk))
if(is.null(cuts)) {
if(length(table(risk))<ngroups) ngroups=length(table(risk))
risk.cat=as.numeric(Hmisc::cut2(risk,g=ngroups))
} else {
risk.cat=as.numeric(Hmisc::cut2(risk,cuts=cuts))
print(table(risk.cat))
ngroups=length(table(risk.cat))
}
emp.risk=numeric(ngroups)
for(i in 1:ngroups){
emp.risk[i]=sum(binary.outcome[risk.cat==i])/sum(weights[risk.cat==i])
}
print(emp.risk)
xx=cbind(risk, "emp.risk"=emp.risk[risk.cat])
xx=xx[order(xx[,"risk"]),]
if (is.null(ylim)) {
if(show.emp.risk) all.=c(emp.risk,risk) else all.=risk
if(scale=="logit") all.=setdiff(all.,0)
ylim=range(all.,na.rm=T)
}
ylab="risk"
if (scale=="logit") {
ylim=logit(ylim)
xx=logit(xx)
ylab="log odd"
}
plot(xx[,"risk"], type="l", xlab="", xaxt="n", main=main, col=1, ylim=ylim, ylab=ylab)
if(show.emp.risk) lines(xx[,"emp.risk"], col=lcol, lty=2)
print(table(xx[,"emp.risk"], useNA="ifany"))
mylegend(col=1:lcol, lty=1:2, x=1, legend=c("model fit","empirical"))
}
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.