Codebase list hmisc / f148276
Update upstream source from tag 'upstream/4.4-2' Update to upstream version '4.4-2' with Debian dir 91ae7b43dac6d8901b9ada4f3a3a4c646096a9b0 Dirk Eddelbuettel 3 years ago
20 changed file(s) with 839 addition(s) and 162 deletion(s). Raw diff Collapse all Expand all
00 Package: Hmisc
1 Version: 4.4-1
2 Date: 2020-08-10
1 Version: 4.4-2
2 Date: 2020-11-25
33 Title: Harrell Miscellaneous
44 Author: Frank E Harrell Jr <fh@fharrell.com>, with
55 contributions from Charles Dupont and many others.
2222 Encoding: UTF-8
2323 RoxygenNote: 7.1.1
2424 NeedsCompilation: yes
25 Packaged: 2020-08-10 12:34:23 UTC; harrelfe
25 Packaged: 2020-11-27 22:05:50 UTC; harrelfe
2626 Repository: CRAN
27 Date/Publication: 2020-08-10 22:02:11 UTC
27 Date/Publication: 2020-11-29 06:00:16 UTC
+19
-14
MD5 less more
00 16f8ea43efba90112351c5ecb33021c4 *COPYING
1 cef04548862dc17bed4924432f389f75 *DESCRIPTION
2 ac280b683a20ca6bf4395c2435fc30de *NAMESPACE
3 20fbb1071259a378a8c11b6c4cf62fd1 *NEWS
1 be0222e063fcd3b8dd31e0fca3c984a3 *DESCRIPTION
2 03aeb100e13bd948cd29ce779bb9da64 *NAMESPACE
3 396507751e208423836df055a792b621 *NEWS
44 76f90acbe9bf8ad2eaaa6b8813eb7c82 *R/AFirst.lib.s
55 62ef16fc333551d6c65f15dedfb5a697 *R/Cs.s
66 85683865508739ff7a19b280e562ab67 *R/GiniMd.s
3232 46e1eda21de21ba281f098f1f9e59bcb *R/describe.s
3333 6cb5b3a77860b605f102bb528e84a071 *R/discrete.s
3434 dcc7f27cea0aa7b110b0eb2be3144d8a *R/dotchart3.s
35 ab7c8cc449c83d6612b125a5da844a5b *R/dotchartpl.s
35 266f05a603d2ddd37ad68ac269a8974d *R/dotchartpl.s
3636 919e94c2c87f720601387c5171b57ffc *R/drawPlot.s
3737 b44d1ac43c71a3a4d3f4820354f78598 *R/ecdf.s
3838 c1e489a07ca933fb02e004489fd3ee8e *R/epi.s
4444 c238614fb72943362407d74442bf236a *R/format.pval.s
4545 aace929824aad6ebdfba910950b6cc2b *R/ftu.s
4646 473bbef2ffd42d37d059d9a4c66fe8a4 *R/gbayes.s
47 7d361bc56d0c8ca82e48b90493e1f0e9 *R/gbayesSeqSim.r
4748 f76f66eae7faef0e445143d9e367619d *R/gettext.s
4849 769b009b031f6f656a0e9b748057bfae *R/ggMisc.s
49 dc8fa7871825b028089e45097975a973 *R/ggfreqScatter.r
50 2d9c334a5b7dddb5f8815e7bc2257c6d *R/ggfreqScatter.r
5051 b70800bb3730a5b84500d43e97d951f4 *R/groupn.s
5152 926950a356c2e54ecd8b2317ec08a5d8 *R/hidingTOC.r
5253 29cf6d5ee0c8444cb58ed0afb85ac57b *R/hist.data.frame.s
8485 1ada10fb83652630cf67024c3f1b421c *R/nobsY.s
8586 27019e9e90730ac0b25014b51d68ec66 *R/nstr.s
8687 cf2866606281d01f39e99d3205f50ae3 *R/num.intercepts.s
88 a50c8048a7546e224ca402d362c084cd *R/pairUpDiff.r
8789 215f9a06ffc6df02fb67576ad0da83b9 *R/panel.abwplot.s
8890 41cf42963c4ceb1cb1233cc1e01b4b99 *R/panel.bpplot.s
8991 1630c1be85251dac0e8cd0243bedd201 *R/pc1.s
9092 25fc3a3579da1d81be1209bea9e63d2b *R/plot.describe.s
9193 c746aa382838bbb7063eec0487f2a64d *R/plotlyM.r
9294 62e6fe8152bc03bb23f83ff687ae4cc2 *R/plsmo.s
93 6c7831f266dc8e510adef7de3be688ef *R/popower.s
95 8393f1ae893c0a37c8d367c1e8cd1c98 *R/popower.s
9496 9a1119c7949407c39968645e921a87a6 *R/print.char.list.s
9597 efe9a9a6c90594d1e7425cc46f2fc061 *R/pstamp.s
9698 620401a9617131be8b2d162fd0466ff7 *R/rcorr.cens.s
97 06948100269e0b2705e6f0b0b1a8d4b9 *R/rcorr.s
99 b9759f7459856f6eeb85e4ee3832f7ad *R/rcorr.s
98100 0c9a682b270e97639d2f7d06c18aaa7f *R/rcorrp.cens.s
99101 62264d939e616b9d1dc7070096b679e5 *R/rcspline.eval.s
100102 bd5869619261cea0beb0bbc69f9d539f *R/rcspline.plot.s
107109 3978f2ee3511549cb906af4618024a53 *R/rm.boot.s
108110 a1763e6429f086dd1783694da9b85342 *R/samplesize.bin.s
109111 130de98e37b80bad7896e6b083be8964 *R/sas.get.s
110 e0871a4ba840ef9892e26282ad7e9bd5 *R/scat1d.s
112 bec8bca50b20965e503b4045b1cd1027 *R/scat1d.s
111113 bcd4fb26ff4be0623165a02d1f874733 *R/score.binary.s
112114 c1639df1f0faf833c9cb440e8367aa53 *R/sedit.s
113115 f361dab6e800a20722ecc51bc08e056d *R/show.pch.s
127129 b54b4f47cd9f797bb3bd49378f0367af *R/summaryM.s
128130 f4e9b0484893f71186d0e527cfcab4b8 *R/summaryP.s
129131 dc27150ba262b8c3844ce46d36f1e492 *R/summaryRc.s
130 f9c4fedbc00b494b1a4322d74e175ee5 *R/summaryS.s
132 0ce190af103865fe51ba98c7bc4dab9d *R/summaryS.s
131133 0af65ac91b7c04b13aac7791bdf758af *R/symbol.freq.s
132134 474107a89a3a998eb78bbee1ac2dfdbe *R/sys.s
133135 0917e77a727d026b5b12511b7710b283 *R/t.test.cluster.s
258260 4c2942aab548faa49c5a389117b48b70 *man/discrete.Rd
259261 94748a0c9bb7d0290a8b2dc7f4a0bb7d *man/dotchart2.Rd
260262 67f28b3925d62e99e612e36532d3149f *man/dotchart3.Rd
261 a716e5e2bf83ff413350835c7058923c *man/dotchartpl.Rd
263 2240c79711c6ab7bc4d204483d3d8258 *man/dotchartpl.Rd
262264 553f6215baa4cad0ac3da850a3cbc2d2 *man/epi.Rd
263265 36c0c98e4e6fcc78979da741c8dfb386 *man/equalBins.Rd
264266 04771822c60d38658c3720cd91bc2113 *man/errbar.Rd
265267 3bc08c67bb196770cb9ff90bdfa15672 *man/escapeRegex.Rd
268 7529c6702e5cd078996c56b78b826c78 *man/estSeqSim.Rd
266269 b8c947ebb9ea63a9ecf543e9248a512b *man/event.chart.Rd
267270 6f9b356b12e05d93c091e033a7d3f4a8 *man/event.convert.Rd
268271 1beee8ab5e6bde656f0daaf7dfe1cbc5 *man/event.history.Rd
270273 6a5adf7af2c0cc119bbcded8040bd34f *man/first.word.Rd
271274 262cfe8e5568594c1ebf3915f81f693a *man/format.df.Rd
272275 041f76e1a74ad1e2361aab071b496dac *man/format.pval.Rd
273 4f6bea3b84c269f61ef89eaf205d8f00 *man/gbayes.Rd
276 ebdd21800ffc4beea994eaac7d5b4668 *man/gbayes.Rd
277 5522826f50906b5ec42f9679d29a5af1 *man/gbayesSeqSim.Rd
274278 65798e0cb002195035db2991427ce63c *man/getHdata.Rd
275279 993b4834e5bbff1af338a24d08cbb6cf *man/getRs.Rd
276280 fa2378b377a46a23ee77876829ebf8c4 *man/getZip.Rd
277281 ac60ee276185df5b27c0819a3e781bae *man/ggMisc.Rd
278 80e90e865489fc64f8414fa058030112 *man/ggfreqScatter.Rd
282 10239b9be9e48b4526e454000124280c *man/ggfreqScatter.Rd
279283 bd30d6597a6591b2d3b6202931502245 *man/hdquantile.Rd
280284 9bde85cef5a8815e6a4adbc2163aa870 *man/hidingTOC.Rd
281285 8256dd2a02e789ad86ec974fa90420cd *man/hist.data.frame.Rd
309313 be66aaec2a64401335915c4671d50918 *man/nobsY.Rd
310314 2f7e82883fd0be5a86a3bec4d4dc250a *man/nstr.Rd
311315 1d560c844c38c8f7289bb27e21efe171 *man/num.intercepts.Rd
316 76b1975307ab89cc2eb19f12040342ce *man/pairUpDiff.Rd
312317 ab8edbb03663977e40ac2ffa3025d39f *man/panel.bpplot.Rd
313318 b8247448b0c2c6d876ce693f23c9fe85 *man/partition.Rd
314319 b34470017ff70767108ad8ee9e169f02 *man/pc1.Rd
315320 17652e750bd0249ced3eac8e145d070d *man/plotCorrPrecision.Rd
316321 88d93c5117d432cbe4ba996c5c74d0d8 *man/plotlyM.Rd
317322 81937c2b724d0ec42f817ebb88748e70 *man/plsmo.Rd
318 53e306cc4774582993e485b474988555 *man/popower.Rd
323 0a2e6d1c460ef2be08a7162c85709e43 *man/popower.Rd
319324 52be83076b1f0442cb5ab5c15bb16a3f *man/print.char.list.Rd
320325 0a25e1c6f349789159b7239ae0f0f620 *man/print.char.matrix.Rd
321326 057da98d466a023bda3991abbf880855 *man/prnz.Rd
334339 0363a389eb98c4c065a4cecae63c7e56 *man/rm.boot.Rd
335340 326928a75e512ba40912a48fc355934e *man/samplesize.bin.Rd
336341 89c4eefbaac065e35e5becf8f6e5319d *man/sasxport.get.Rd
337 c7247e422a999499fc292af52da803a8 *man/scat1d.Rd
342 ca36423dc97de82f1aef3c6d260c7dd8 *man/scat1d.Rd
338343 db05fb7a0b4ab6f2e3296104ab70dca1 *man/score.binary.Rd
339344 08744258b9fd7ca14e1b32ba62ca245f *man/sedit.Rd
340345 8dbcd26d8e34394832eaa6c2f32931a3 *man/show.pch.Rd
0 export("%nin%",abs.error.pred,addMarginal,all.digits,all.is.numeric,approxExtrap,areg,areg.boot,aregImpute,aregTran,arrGrob,as.discrete,asNumericMatrix,ballocation,bezier,binconf,biVar,bootkm,bpower,bpower.sim,bpplot,bpplotM,bpplt,bppltp,bpx,bsamsize,bystats,bystats2,capitalize,catTestchisq,Cbind,ceil,character.table,chiSquare,ciapower,cleanup.import,clowess,cnvrt.coords,code.levels,colorFacet,combine.levels,combineLabels,combplotp,confbar,consolidate,contents,conTestkw,cpower,Cs,csv.get,cumcategory,curveRep,curveSmooth,cut2,datadensity,dataDensityString,dataframeReduce,dataRep,ddmmmyy,deff,describe,discrete,dhistboxp,dotchart2,dotchart3,dotchartp,dotchartpl,Dotplot,drawPlot,dvi,dvigv,dvips,ecdfpM,Ecdf,equalBins,errbar,escapeBS,escapeRegex,event.chart,event.convert,event.history,expr.tree,fillin,find.matches,first.word,fit.mult.impute,format.df,format.pval,formatCats,formatCons,formatDateTime,formatdescribeSingle,formatSep,formatTestStats,ftupwr,ftuss,Function,gbayes,gbayes1PowerNP,gbayes2,gbayesMixPost,gbayesMixPowerNP,gbayesMixPredNoData,get2rowHeads,getHdata,getLatestSource,GetModelFrame,getRs,getZip,ggfreqScatter,GiniMd,Gompertz2,groupn,grType,hdquantile,hidingTOC,hist.data.frame,histbackback,histboxp,histboxpM,histSpike,histSpikeg,hoeffd,html,htmlGreek,htmlSN,htmlSpecial,htmlSpecialType,htmlTranslate,htmlVerbatim,importConvertDateTime,improveProb,impute,impute.transcan,inmChoice,inverseFunction,invertTabulated,is.discrete,is.imputed,is.mChoice,is.present,is.special.miss,james.stein,jitter2,keepHattrib,Key,Key2,km.quick,knitrSet,labcurve,label,Label,labelLatex,labelPlotmath,Lag,largest.empty,latex,latexBuild,latexCheckOptions,latexDotchart,latexNeedle,latexSN,latexTabular,latexTherm,latexTranslate,latexVerbatim,list.tree,llist,lm.fit.qr.bare,Load,Lognorm2,logrank,lookupSASContents,lrcum,makeNames,makeNstr,makeSteps,mApply,markupSpecs,mask,match.mChoice,matchCases,matrix2dataFrame,matxv,mbarclPanel,mbarclpl,mChoice,mdb.get,Mean,medvPanel,medvpl,Merge,mgp.axis,mgp.axis.labels,mhgr,minor.tick,monotone,monthDays,mtitle,multEventChart,multLines,na.delete,na.detail.response,na.include,na.keep,na.pattern,na.retain,naclus,nafitted.delete,Names2names,naplot,napredict.delete,naprint.delete,naprint.keep,naresid.delete,naresid.keep,nFm,nobsY,nomiss,nstr,num.denom.setup,num.intercepts,numeric.string,numericScale,oPar,optionsCmds,ordGridFun,ordTestpo,outerText,panel.bpplot,panel.Dotplot,panel.Ecdf,panel.plsmo,panel.xYplot,parGrid,partition.matrix,partition.vector,pasteFit,pBlock,pc1,plotCorrPrecision,plotp,plotlyM,plotlyParm,plotlySave,plotmathTranslate,plotMultSim,plotpsummaryM,plsmo,pngNeedle,pomodm,popower,posamsize,prepanel.Dotplot,prepanel.Ecdf,prepanel.xYplot,print.char.matrix,prList,prn,propsPO,propsTrans,prselect,prType,pstamp,putHcap,putHfig,putKey,putKeyEmpty,Quantile,Quantile2,rcorr,rcorr.cens,rcorrcens,rcorrp.cens,rcspline.eval,rcspline.plot,rcspline.restate,rcsplineFunction,read.xportDataload,readSAScsv,redun,reformM,reLabelled,replace.substring.wild,reShape,responseSummary,restoreHattrib,rlegend,rlegendg,rm.boot,rMultinom,roundN,roundPOSIXt,samplesize.bin,sas.codes,sas.get,sas.get.macro,sasdsLabels,sasxport.get,Save,scat1d,score.binary,sedit,sepUnitsTrans,setParNro,setTrellis,show.col,show.pch,showPsfrag,simplifyDims,simRegOrd,sKey,smean.cl.boot,smean.cl.normal,smean.sd,smean.sdl,smearingEst,smedian.hilow,solvet,somers2,spearman,spearman.test,spearman2,spower,spss.get,src,stat_plsmo,stata.get,StatPlsmo,stepfun.eval,stratify,strgraphwrap,string.bounding.box,string.break.line,stringDims,stripChart,subplot,substi,substi.source,substring.location,substring2,summarize,summaryD,summaryDp,summaryM,summaryP,summaryRc,summaryS,symbol.freq,sys,t.test.cluster,table_formatpct,table_freq,table_latexdefs,table_N,table_pc,table_trio,tabulr,termsDrop,testDateTime,tex,tobase64image,transace,transcan,translate,trap.rule,trellis.strip.blank,truncPOSIXt,uncbind,units,unPaste,upData,upFirst,valueLabel,valueName,valueTags,valueUnit,var.inner,varclus,Weibull2,whichClosek,whichClosePW,whichClosest,wtd.Ecdf,wtd.loess.noiter,wtd.mean,wtd.quantile,wtd.rank,wtd.table,wtd.var,xInch,xless,xy.group,xYplot,xySortNoDupNoNA,yearDays,yInch,ynbind,zoom,"[<-.discrete","consolidate<-","is.na<-.discrete","label<-","length<-.discrete","substring2<-","valueLabel<-","valueName<-","valueTags<-","valueUnit<-")
0 export("%nin%",abs.error.pred,addMarginal,all.digits,all.is.numeric,approxExtrap,areg,areg.boot,aregImpute,aregTran,arrGrob,as.discrete,asNumericMatrix,ballocation,bezier,binconf,biVar,bootkm,bpower,bpower.sim,bpplot,bpplotM,bpplt,bppltp,bpx,bsamsize,bystats,bystats2,capitalize,catTestchisq,Cbind,ceil,character.table,chiSquare,ciapower,cleanup.import,clowess,cnvrt.coords,code.levels,colorFacet,combine.levels,combineLabels,combplotp,confbar,consolidate,contents,conTestkw,cpower,Cs,csv.get,cumcategory,curveRep,curveSmooth,cut2,datadensity,dataDensityString,dataframeReduce,dataRep,ddmmmyy,deff,describe,discrete,dhistboxp,dotchart2,dotchart3,dotchartp,dotchartpl,Dotplot,drawPlot,dvi,dvigv,dvips,ecdfpM,Ecdf,equalBins,errbar,escapeBS,escapeRegex,estSeqSim,event.chart,event.convert,event.history,expr.tree,fillin,find.matches,first.word,fit.mult.impute,format.df,format.pval,formatCats,formatCons,formatDateTime,formatdescribeSingle,formatSep,formatTestStats,ftupwr,ftuss,Function,gbayes,gbayes1PowerNP,gbayes2,gbayesMixPost,gbayesMixPowerNP,gbayesMixPredNoData,gbayesSeqSim,get2rowHeads,getHdata,getLatestSource,GetModelFrame,getRs,getZip,ggfreqScatter,GiniMd,Gompertz2,groupn,grType,hdquantile,hidingTOC,hist.data.frame,histbackback,histboxp,histboxpM,histSpike,histSpikeg,hoeffd,html,htmlGreek,htmlSN,htmlSpecial,htmlSpecialType,htmlTranslate,htmlVerbatim,importConvertDateTime,improveProb,impute,impute.transcan,inmChoice,inverseFunction,invertTabulated,is.discrete,is.imputed,is.mChoice,is.present,is.special.miss,james.stein,jitter2,keepHattrib,Key,Key2,km.quick,knitrSet,labcurve,label,Label,labelLatex,labelPlotmath,Lag,largest.empty,latex,latexBuild,latexCheckOptions,latexDotchart,latexNeedle,latexSN,latexTabular,latexTherm,latexTranslate,latexVerbatim,list.tree,llist,lm.fit.qr.bare,Load,Lognorm2,logrank,lookupSASContents,lrcum,makeNames,makeNstr,makeSteps,mApply,markupSpecs,mask,match.mChoice,matchCases,matrix2dataFrame,matxv,mbarclPanel,mbarclpl,mChoice,mdb.get,Mean,medvPanel,medvpl,Merge,mgp.axis,mgp.axis.labels,mhgr,minor.tick,monotone,monthDays,mtitle,multEventChart,multLines,na.delete,na.detail.response,na.include,na.keep,na.pattern,na.retain,naclus,nafitted.delete,Names2names,naplot,napredict.delete,naprint.delete,naprint.keep,naresid.delete,naresid.keep,nFm,nobsY,nomiss,nstr,num.denom.setup,num.intercepts,numeric.string,numericScale,oPar,optionsCmds,ordGridFun,ordTestpo,outerText,pairUpDiff,panel.bpplot,panel.Dotplot,panel.Ecdf,panel.plsmo,panel.xYplot,parGrid,partition.matrix,partition.vector,pasteFit,pBlock,pc1,plotCorrPrecision,plotp,plotlyM,plotlyParm,plotlySave,plotmathTranslate,plotMultSim,plotpsummaryM,plsmo,pngNeedle,pomodm,popower,posamsize,prepanel.Dotplot,prepanel.Ecdf,prepanel.xYplot,print.char.matrix,prList,prn,propsPO,propsTrans,prselect,prType,pstamp,putHcap,putHfig,putKey,putKeyEmpty,Quantile,Quantile2,rcorr,rcorr.cens,rcorrcens,rcorrp.cens,rcspline.eval,rcspline.plot,rcspline.restate,rcsplineFunction,read.xportDataload,readSAScsv,redun,reformM,reLabelled,replace.substring.wild,reShape,responseSummary,restoreHattrib,rlegend,rlegendg,rm.boot,rMultinom,roundN,roundPOSIXt,samplesize.bin,sas.codes,sas.get,sas.get.macro,sasdsLabels,sasxport.get,Save,scat1d,score.binary,sedit,sepUnitsTrans,setParNro,setTrellis,show.col,show.pch,showPsfrag,simplifyDims,simRegOrd,sKey,smean.cl.boot,smean.cl.normal,smean.sd,smean.sdl,smearingEst,smedian.hilow,solvet,somers2,spearman,spearman.test,spearman2,spower,spss.get,src,stat_plsmo,stata.get,StatPlsmo,stepfun.eval,stratify,strgraphwrap,string.bounding.box,string.break.line,stringDims,stripChart,subplot,substi,substi.source,substring.location,substring2,summarize,summaryD,summaryDp,summaryM,summaryP,summaryRc,summaryS,symbol.freq,sys,t.test.cluster,table_formatpct,table_freq,table_latexdefs,table_N,table_pc,table_trio,tabulr,termsDrop,testDateTime,tex,tobase64image,transace,transcan,translate,trap.rule,trellis.strip.blank,truncPOSIXt,uncbind,units,unPaste,upData,upFirst,valueLabel,valueName,valueTags,valueUnit,var.inner,varclus,Weibull2,whichClosek,whichClosePW,whichClosest,wtd.Ecdf,wtd.loess.noiter,wtd.mean,wtd.quantile,wtd.rank,wtd.table,wtd.var,xInch,xless,xy.group,xYplot,xySortNoDupNoNA,yearDays,yInch,ynbind,zoom,"[<-.discrete","consolidate<-","is.na<-.discrete","label<-","length<-.discrete","substring2<-","valueLabel<-","valueName<-","valueTags<-","valueUnit<-")
11
22 useDynLib(Hmisc, .registration=TRUE, .fixes="F_")
33
0 Changes in version 4.4-2 (2020-11-25)
1 * rcorr: captured loss of precision leading to square root of a negative number. Thanks: Ann Voss <avoss@novanthealth.org>
2 * summaryS: sapply was doubling names
3 * pairUpDiff: created for dotchartpl - function to pair up grouped observations for sorting by descending differences and computing approximate confidence intervals for the difference given individual confidence limits
4 * dotchartpl: added new arguments and functionality to handle quantities other than proportions (e.g., hazards) and streamlined code using pairUpDiff
5 * propsPO: added tooltip text for gpplot that will be transmitted to ggplotly; reversed order of categories so lowest category put at bottom of bar
6 * dotchartpl: suppressed confidence intervals when the number of groups is not 2; fixed bug where hover text confidence intervals were repeats of the last computed interval instead of properly using all the intervals; added dec argument
7 * added estSeqSim and gbayesSeqSim functions
8 * ggfreqScatter: stopped varying alpha and just varied color, as alpha collided with the color scheme
9 * histSpike: added minimal=TRUE and bins arguments
10
011 Changes in version 4.4-1 (2020-08-07)
112 * popower: added approximate S.E. of log(OR) in results
213 * propsPO: new function for exploring proportional odds
0 dotchartpl <- function(x, major=NULL, minor=NULL, group=NULL, mult=NULL,
0 dotchartpl <- function(x, major=NULL, minor=NULL,
1 group=NULL, mult=NULL,
12 big=NULL, htext=NULL,
23 num=NULL, denom=NULL,
4 numlabel='', denomlabel='',
5 fun=function(x) x, ifun=function(x) x,
6 op='-',
37 lower=NULL, upper=NULL,
48 refgroup=NULL, sortdiff=TRUE, conf.int=0.95,
59 minkeep=NULL,
610 xlim=NULL, xlab='Proportion',
711 tracename=NULL, limitstracename='Limits',
812 nonbigtracename='Stratified Estimates',
9 width=800, height=NULL,
13 dec=3, width=800, height=NULL,
1014 col=colorspace::rainbow_hcl
1115 ) {
16
1217 mu <- markupSpecs$html
1318 bold <- mu$bold
1419
1621 max(c(x, upper), na.rm=TRUE))
1722
1823 majorpres <- length(major) > 0
19 major <- if(majorpres) as.character(major) else rep('', length(x))
24 major <- if(majorpres) as.character(major) else rep(' ', length(x))
2025 minorpres <- length(minor) > 0
26 if(! (majorpres || minorpres)) stop('must specify major or minor or both')
27
2128 grouppres <- length(group) > 0 ## superpositioning variable
2229 multpres <- length(mult) > 0
2330 limspres <- length(lower) * length(upper) > 0
2431 rgpres <- length(refgroup) > 0
2532
26 if(minorpres) minor <- as.character(minor)
27 if(grouppres) group <- as.character(group)
28 if(multpres) mult <- as.character(mult)
29 ugroup <- if(grouppres) unique(group)
33 if(minorpres) minor <- as.character(minor)
34 if(grouppres) group <- as.character(group)
35 if(multpres) mult <- as.character(mult)
36 ugroup <- if(grouppres) unique(group) else ''
37
38 fmt <- function(x) format(round(x, dec))
3039
3140 if(rgpres) {
3241 if(! grouppres || multpres || length(big))
4857 col(1) }
4958 if(! length(col) && ! grouppres) cols <- 'black'
5059
51 levelsRemoved <- character(0)
52
53 if(rgpres && (sortdiff || length(minkeep)) && minorpres) {
54 ## For each major x minor subgroup, hold on to x for both groups
55
56 g <- function(i, wcat) {
57 gr <- group[i]
58 if(any(gr == wcat)) x[i[gr == wcat]] else NA
59 }
60
61 xa <- tapply(1 : length(x), llist(major, minor), g, wcat=refgroup)
62 xb <- tapply(1 : length(x), llist(major, minor), g, wcat=altgroup)
63
64 ## Bug in x[cbind()] if major is constant
65 if(! majorpres) {
66 xa <- xa[, minor]
67 xb <- xb[, minor]
68 }
69 else {
70 xa <- xa[cbind(major, minor)]
71 xb <- xb[cbind(major, minor)]
72 }
73
74 if(length(minkeep)) {
75 w <- if(majorpres) paste0(major, ':', minor) else minor
76 levelsRemoved <- unique(w[! is.na(xa + xb) &
77 xa < minkeep & xb < minkeep])
78 j <- which(! is.na(xa + xb) & (xa >= minkeep | xb >= minkeep))
79 }
80 else
81 j <- 1 : length(x)
82
83 ## Sort minor categories by descending order of between-group differences
84
85 i <- order(major, ifelse(is.na(xa + xb), -Inf, - (xb - xa)))
86 i <- intersect(i, j)
60 dropped <- character(0)
61 D <- NULL
62
63 if(rgpres && minorpres) {
64
65 z <- pairUpDiff(x, major, minor, group,
66 refgroup=refgroup, lower=lower, upper=upper,
67 minkeep=minkeep, sortdiff=sortdiff, conf.int=conf.int)
68
69 Z <- z$X
70 D <- z$D
71 dropped <- z$dropped
72 i <- Z[,'subscripts']
73
8774 x <- x[i]
8875 major <- major[i]
8976 minor <- minor[i]
9582 if(length(lower)) lower <- lower[i]
9683 if(length(upper)) upper <- upper[i]
9784 if(length(htext)) htext <- htext[i]
98 }
85 } # end if(rgpres && ...
9986
10087
10188 ht <- htext
89 if(numlabel != '') numlabel <- paste0(' ', numlabel)
90 if(denomlabel != '') denomlabel <- paste0(' ', denomlabel)
10291 if(length(num))
10392 ht <- paste0(ht, if(length(htext)) mu$lspace,
104 round(x, 3), mu$lspace,
105 mu$frac(num, denom, size=95))
93 fmt(fun(x)), mu$lspace,
94 mu$frac(paste0(num, numlabel),
95 paste0(denom, denomlabel), size=95))
96
97 ## if confidence interval for differences are not to be displayed,
98 ## put point estimate confidence intervals in point hypertext
99 if(length(ugroup) != 2 && limspres)
100 ht <- paste0(ht, ' [',
101 fmt(fun(lower)), ', ',
102 fmt(fun(upper)), ']')
103
106104 ht <- paste0(ht, if(length(ht)) '<br>',
107105 if(majorpres) paste0(major, ': '))
108106 if(minorpres) ht <- paste0(ht, minor)
148146 y <- y - 0.4
149147 lines <- lines + (mi != '')
150148 j <- which(major == ma & minor == mi)
151 m <- length(j)
152149 X <- c(X, x[j])
153150 Y <- c(Y, ifelse(big[j], y, y - .14))
154 if(rgpres) {
155 ja <- which(major == ma & minor == mi & group == refgroup)
156 jb <- which(major == ma & minor == mi & group == altgroup)
157 diff <- x[jb] - x[ja]
151 if(length(D)) {
152 k <- which(D$major == ma & D$minor == mi)
153 if(length(k) != 1)
154 stop('must have one observation for a major/minor/group combination')
155 diff <- D$diff[k]
158156 coldiff <- c(coldiff,
159157 ifelse(diff > 0,
160158 paste0(altgroup, ' ',
165163
166164 htd <- if(majorpres) paste0(ma, ': ') else ''
167165 if(minorpres) htd <- paste0(htd, mi)
168
169 htd <- paste0(htd, '<br>', altgroup, ' - ', refgroup, ': ',
170 round(diff, 3))
171
172 if(! is.logical(conf.int) && length(num) && length(denom)) {
173 zcrit <- qnorm((1 + conf.int) / 2)
174 ## If either proportion is 0 or 1, backsolve variance from
175 ## Wilson interval
176 wi <- function(p, n) {
177 if(p == 0 || p == 1) {
178 ci <- binconf(round(n * p), n)[1, ]
179 w <- ci['Upper'] - ci['Lower']
180 (0.5 * w / zcrit) ^ 2
181 } else p * (1. - p) / n
182 }
183 va <- wi(x[ja], denom[ja])
184 vb <- wi(x[jb], denom[jb])
185 se <- sqrt(va + vb)
186 dlower <- format(round(x[jb] - x[ja] - zcrit * se, 3), nsmall=3)
187 dupper <- format(round(x[jb] - x[ja] + zcrit * se, 3), nsmall=3)
166
167 htd <- paste0(htd, '<br>', altgroup, ' ', op, ' ',
168 refgroup, ': ', fmt(fun(diff)))
169 if(! is.logical(conf.int) && limspres && length(D)) {
170 dlower <- fmt(fun(D$lower[k]))
171 dupper <- fmt(fun(D$upper[k]))
188172 htd <- paste0(htd, '<br>', conf.int, ' C.L.: [', dlower,
189173 ', ', dupper, ']')
190 Diff <- c(Diff, x[jb] - x[ja])
191 difflower <- c(difflower, (x[jb] + x[ja]) / 2 - 0.5 * zcrit * se)
192 diffupper <- c(diffupper, (x[jb] + x[ja]) / 2 + 0.5 * zcrit * se)
174 Diff <- c(Diff, diff)
175 difflower <- c(difflower, D$lowermid[k])
176 diffupper <- c(diffupper, D$uppermid[k])
193177 }
194178 htdiff <- c(htdiff, htd)
195 }
196 if(limspres) {
179 } # end if(length(D))
180 if(limspres && ! length(D)) {
197181 Lower <- c(Lower, lower[j])
198182 Upper <- c(Upper, upper[j])
199 Htextl <- paste0('[', format(lower[j], digits=4), ', ',
200 format(upper[j], digits=4), ']')
183 Htextl <- c(Htextl,
184 paste0('[',
185 fmt(fun(lower[j])), ', ',
186 fmt(fun(upper[j])), ']' ) )
201187 }
202188 Group <- c(Group, group[j])
203 Big <- c(Big, big[j])
189 Big <- c(Big, big[j])
204190 Htext <- c(Htext, ht[j])
205191
206192 if(mi != '') {
213199
214200 d <- data.frame(X, Y, Group, Htext, Big)
215201
216 if(limspres) {
202 if(limspres && ! length(D)) {
217203 d$Lower <- Lower
218204 d$Upper <- Upper
219205 d$Htextl <- Htextl
220206 }
221
222 if(! grouppres) d$Group <- NULL ####
223 if(any(d$Big)) {
207
208 if(! grouppres) d$Group <- NULL
209
210 if(any(d$Big)) {
224211 db <- subset(d, Big) # non-stratified estimates
225212 ## For some reason, colors= in add_marker did not always take
226213 p <- plotly::plot_ly(data=db, colors=cols,
227214 height=if(length(height)) height else
228215 plotlyParm$heightDotchart(lines), width=width)
229
230 if(length(difflower)) {
216
217 if(limspres && length(D)) {
231218 ddiff <- data.frame(Diff, difflower, diffupper, Ydiff, coldiff, htdiff)
219
232220 ## Could not get color=coldiff to work; perhaps conflicted with
233221 ## earlier use of color =
222
234223 if(any(Diff > 0))
235224 p <- plotly::add_segments(p, data=subset(ddiff, Diff > 0),
236225 x= ~ difflower, xend= ~ diffupper,
240229 name = paste0(htmlSpecial('half'),
241230 ' CL of difference<br>',
242231 coldiff[Diff > 0][1]))
232
243233 if(any(Diff <= 0))
244234 p <- plotly::add_segments(p, data=subset(ddiff, Diff <= 0),
245235 x= ~ difflower, xend= ~ diffupper,
250240 ' CL of difference<br>',
251241 coldiff[Diff <= 0][1]))
252242 }
243
253244 ## tracename and limitstracename are used if groups not used
254 if(limspres)
245
246 if(limspres && ! length(D) && length(ugroup) == 2)
255247 p <- if(grouppres)
256 plotly::add_segments(p, x=~ Lower, xend=~ Upper,
248 plotly::add_segments(p, data=db,
249 x=~ Lower, xend=~ Upper,
257250 y=~ Y, yend=~ Y,
258251 color=~ Group, text= ~ Htextl,
259252 colors=cols, hoverinfo='text')
263256 text= ~ Htextl,
264257 color=I('lightgray'), hoverinfo='text',
265258 name=limitstracename)
266
259
267260 p <- if(grouppres)
268261 plotly::add_markers(p, x=~ X, y=~ Y,
269262 color=~ Group, text=~ Htext,
275268 name=if(length(tracename)) tracename
276269 else
277270 if(any(! d$Big)) 'All' else '')
278
271
279272 } else p <- plotly::plot_ly(colors=cols) # Big is not used
280273
281274 if(any(! d$Big)) {
282275 dnb <- subset(d, ! Big) # stratified estimates
283276 if(limspres)
284 p <- if(grouppres)
277 p <- if(grouppres && length(ugroup) == 2)
285278 plotly::add_segments(p, data=dnb,
286279 x=~ Lower, xend=~ Upper,
287280 y=~ Y, yend=~ Y,
314307 marker=list(opacity=0.45, size=4),
315308 color=I('black'), hoverinfo='text',
316309 name=nonbigtracename)
317
318 } # end if(any(! Big))
310
311 } # end if(any(! Big))
312
319313 leftmargin <- plotlyParm$lrmargin(ytnb)
314 tickvals <- pretty(fun(xlim), n=10)
315 xaxis <- list(title=xlab, range=xlim, zeroline=FALSE,
316 tickvals=ifun(tickvals), ticktext=tickvals)
317 yaxis <- list(title='',
318 range=c(min(Y) - 0.2, 0.2),
319 zeroline=FALSE, tickvals=yl, ticktext=yt)
320
320321 p <- plotly::layout(p,
321 xaxis=list(title=xlab,
322 range=xlim,
323 zeroline=FALSE),
324 yaxis=list(title='',
325 range=c(min(Y) - 0.2, 0.2),
326 zeroline=FALSE, tickvals=yl, ticktext=yt),
322 xaxis=xaxis,
323 yaxis=yaxis,
327324 margin=list(l=leftmargin),
328325 legend=list(traceorder=if(length(difflower))
329326 'reversed' else 'normal'))
330327
331 attr(p, 'levelsRemoved') <- levelsRemoved
328 attr(p, 'levelsRemoved') <- dropped
332329 p
333330 }
0 ##' Simulate Bayesian Sequential Treatment Comparisons Using a Gaussian Model
1 ##'
2 ##' Simulate a sequential trial under a Gaussian model for parameter estimates, and Gaussian priors using simulated estimates and variances returned by `estSeqSim`. For each row of the data frame `est` and for each prior/assertion combination, computes the posterior probability of the assertion.
3 ##' @title gbayesSeqSim
4 ##' @param est data frame created by `estSeqSim()`
5 ##' @param asserts list of lists. The first element of each list is the user-specified name for each assertion/prior combination, e.g., `"efficacy"`. The other elements are, in order, a character string equal to "<", ">", or "in", a parameter value `cutoff` (for "<" and ">") or a 2-vector specifying an interval for "in", and either a prior distribution mean and standard deviation named `mu` and `sigma` respectively, or a parameter value (`"cutprior"`) and tail area `"tailprob"`. If the latter is used, `mu` is assumed to be zero and `sigma` is solved for such that P(parameter > 'cutprior') = P(parameter < - 'cutprior') = `tailprob`.
6 ##' @return a data frame with number of rows equal to that of `est` with a number of new columns equal to the number of assertions added. The new columns are named `p1`, `p2`, `p3`, ... (posterior probabilities), `mean1`, `mean2`, ... (posterior means), and `sd1`, `sd2`, ... (posterior standard deviations). The returned data frame also has an attribute `asserts` added which is the original `asserts` augmented with any derived `mu` and `sigma` and converted to a data frame, and another attribute `alabels` which is a named vector used to map `p1`, `p2`, ... to the user-provided labels in `asserts`.
7 ##' @author Frank Harrell
8 ##' @seealso `gbayes()`, `estSeqSim()`
9 ##' @examples
10 ##' \dontrun{
11 ##' # Simulate Bayesian operating characteristics for an unadjusted
12 ##' # proportional odds comparison (Wilcoxon test)
13 ##' # For 100 simulations, 5 looks, 2 true parameter values, and
14 ##' # 2 assertion/prior combinations, compute the posterior probability
15 ##' # Use a low-level logistic regression call to speed up simuluations
16 ##' # Use data.table to compute various summary measures
17 ##' # Total simulation time: 2s
18 ##' lfit <- function(x, y) {
19 ##' f <- rms::lrm.fit(x, y)
20 ##' k <- length(coef(f))
21 ##' c(coef(f)[k], vcov(f)[k, k])
22 ##' }
23 ##' gdat <- function(beta, n1, n2) {
24 ##' # Cell probabilities for a 7-category ordinal outcome for the control group
25 ##' p <- c(2, 1, 2, 7, 8, 38, 42) / 100
26 ##'
27 ##' # Compute cell probabilities for the treated group
28 ##' p2 <- pomodm(p=p, odds.ratio=exp(beta))
29 ##' y1 <- sample(1 : 7, n1, p, replace=TRUE)
30 ##' y2 <- sample(1 : 7, n2, p2, replace=TRUE)
31 ##' list(y1=y1, y2=y2)
32 ##' }
33 ##'
34 ##' # Assertion 1: log(OR) < 0 under prior with prior mean 0.1 and sigma 1 on log OR scale
35 ##' # Assertion 2: OR between 0.9 and 1/0.9 with prior mean 0 and sigma computed so that
36 ##' # P(OR > 2) = 0.05
37 ##' asserts <- list(list('Efficacy', '<', 0, mu=0.1, sigma=1),
38 ##' list('Similarity', 'in', log(c(0.9, 1/0.9)),
39 ##' cutprior=log(2), tailprob=0.05))
40 ##'
41 ##' set.seed(1)
42 ##' est <- estSeqSim(c(0, log(0.7)), looks=c(50, 75, 95, 100, 200),
43 ##' gendat=gdat,
44 ##' fitter=lfit, nsim=100)
45 ##' z <- gbayesSeqSim(est, asserts)
46 ##' head(z)
47 ##' attr(z, 'asserts')
48 ##'
49 ##' # Compute the proportion of simulations that hit targets (different target posterior
50 ##' # probabilities for efficacy vs. similarity)
51 ##'
52 ##' # For the efficacy assessment compute the first look at which the target
53 ##' # was hit (set to infinity if never hit)
54 ##' require(data.table)
55 ##' z <- data.table(z)
56 ##' u <- z[, .(first=min(p1 > 0.95])), by=.(parameter, sim)]
57 ##' # Compute the proportion of simulations that ever hit the target and
58 ##' # that hit it by the 100th subject
59 ##' u[, .(ever=mean(first < Inf)), by=.(parameter)]
60 ##' u[, .(by75=mean(first <= 100)), by=.(parameter)]
61 ##' }
62 ##' @md
63 gbayesSeqSim <- function(est, asserts) {
64 nas <- length(asserts)
65 alabels <- character(nas)
66 for(i in 1 : nas) {
67 a <- asserts[[i]]
68 nam <- names(a)
69 if(nam[1] == '') nam[1] <- 'label'
70 if(nam[2] == '') nam[2] <- 'dir'
71 if(nam[3] == '') nam[3] <- 'cutoff'
72 names(a) <- nam
73 if((length(a$cutoff) == 2) != (a$dir == 'in'))
74 stop('mismatch of direction and length of cutoff in asserts')
75 if(any(c('mu', 'sigma') %nin% names(a))) {
76 a$mu <- 0
77 a$sigma <- - a$cutprior / qnorm(a$tailprob)
78 if(a$sigma <= 0) stop('error in specification of cutoff or tailprob')
79 }
80 w <- format(round(a$cutoff, 3))
81 a$assertion <- paste(if(a$dir %in% c('<','>')) a$dir,
82 if(length(a$cutoff) == 2)
83 paste0('[', w[1], ',', w[2], ']')
84 else
85 w[1])
86 asserts[[i]] <- a
87 alabels[i] <- a$label
88 }
89
90 N <- nrow(est)
91
92 ests <- est$est
93 vests <- est$vest
94 ## For each simulated parameter estimate and variance compute nas
95 ## posterior probabilities
96
97 for(i in 1 : nas) {
98 a <- asserts[[i]]
99 dir <- a$dir
100 cutoff <- a$cutoff
101 mu <- a$mu
102 sigma2 <- (a$sigma) ^ 2
103 var.post <- 1. / (1. / sigma2 + 1. / vests)
104 mean.post <- (mu / sigma2 + ests / vests) * var.post
105 sd.post <- sqrt(var.post)
106
107 pp <-
108 switch(dir,
109 '<' = pnorm(cutoff, mean.post, sd.post),
110 '>' = pnorm(cutoff, mean.post, sd.post, lower.tail=FALSE),
111 'in' = pnorm(cutoff[2], mean.post, sd.post) -
112 pnorm(cutoff[1], mean.post, sd.post))
113 label(pp) <- a$label
114 est[[paste0('p', i)]] <- pp
115 est[[paste0('mean', i)]] <- mean.post
116 est[[paste0('sd', i)]] <- sd.post
117 }
118
119 a <- asserts
120 w <- NULL
121 for(j in 1 : nas) {
122 a <- asserts[[j]]
123 a$dir <- a$cutoff <- NULL
124 if(any(c('cutprior', 'tailprob') %nin% names(a)))
125 a$cutprior <- a$tailprob <- NA
126 w <- rbind(w, as.data.frame(a))
127 }
128 attr(est, 'asserts') <- w
129 names(alabels) <- paste0('p', 1 : nas)
130 attr(est, 'alabels') <- alabels
131 est
132 }
133
134 ##' Simulate Comparisons For Use in Sequential Clinical Trial Simulations
135 ##'
136 ##' Simulates sequential clinical trials. Looks are done sequentially at observation numbers given in the vector `looks` with the earliest possible look being at observation 2. For each true effect parameter value, simulation, and at each look, runs a function to compute the estimate of the parameter of interest along with its variance. For each simulation, data are first simulated for the last look, and these data are sequentially revealed for earlier looks. The user provides a function `gendat` that given a true effect of `parameter` and the two sample sizes (for treatment groups 1 and 2) returns a list with vectors `y1` and `y2` containing simulated data. The user also provides a function `fitter` with arguments `x` (group indicator 0/1) and `y` (response variable) that returns a 2-vector containing the effect estimate and its variance. `parameter` is usually on the scale of a regression coefficient, e.g., a log odds ratio.
137 ##' @title estSeqSim
138 ##' @param parameter vector of true parameter (effects; group differences) values
139 ##' @param looks integer vector of observation numbers at which posterior probabilities are computed
140 ##' @param gendat a function of three arguments: true parameter value (scalar), sample size for first group, sample size for second group
141 ##' @param fitter a function of two arguments: 0/1 group indicator vector and the dependent variable vector
142 ##' @param nsim number of simulations (default is 1)
143 ##' @param progress set to `TRUE` to send current iteration number to the console
144 ##' @return a data frame with number of rows equal to the product of `nsim`, the length of `looks`, and the length of `parameter`.
145 ##' @author Frank Harrell
146 ##' @seealso `gbayesSeqSim()`
147 ##' @examples
148 ##' # Run 100 simulations, 5 looks, 2 true parameter values
149 ##' # Total simulation time: 2s
150 ##' lfit <- function(x, y) {
151 ##' f <- rms::lrm.fit(x, y)
152 ##' k <- length(coef(f))
153 ##' c(coef(f)[k], vcov(f)[k, k])
154 ##' }
155 ##' gdat <- function(beta, n1, n2) {
156 ##' # Cell probabilities for a 7-category ordinal outcome for the control group
157 ##' p <- c(2, 1, 2, 7, 8, 38, 42) / 100
158 ##'
159 ##' # Compute cell probabilities for the treated group
160 ##' p2 <- pomodm(p=p, odds.ratio=exp(beta))
161 ##' y1 <- sample(1 : 7, n1, p, replace=TRUE)
162 ##' y2 <- sample(1 : 7, n2, p2, replace=TRUE)
163 ##' list(y1=y1, y2=y2)
164 ##' }
165 ##'
166 ##' set.seed(1)
167 ##' est <- estSeqSim(c(0, log(0.7)), looks=c(50, 75, 95, 100, 200),
168 ##' gendat=gdat,
169 ##' fitter=lfit, nsim=100)
170 ##' head(est)
171 ##' @md
172 estSeqSim <- function(parameter, looks, gendat, fitter, nsim=1,
173 progress=FALSE) {
174
175 looks <- sort(looks)
176 nlook <- length(looks)
177 N <- max(looks)
178 np <- length(parameter)
179 nc <- nsim * nlook * np
180 parm <- est <- vest <- numeric(nc)
181 look <- sim <- integer(nc)
182
183 ## For each simulation and each parameter value, simulate data for the
184 ## whole study
185
186 is <- 0
187 for(isim in 1 : nsim) {
188 if(progress) cat('Simulation', isim, '\r')
189 for(param in parameter) {
190 X <- sample(0 : 1, N, replace=TRUE)
191 dat <- gendat(param, sum(X == 0), sum(X == 1))
192 Y <- rep(NA, N)
193 Y[X == 0] <- dat$y1
194 Y[X == 1] <- dat$y2
195
196 ## For each look compute the parameter estimate and its variance
197
198 for(l in looks) {
199 f <- fitter(X[1 : l], Y[1 : l])
200 is <- is + 1
201 sim[is] <- isim
202 parm[is] <- param
203 look[is] <- l
204 est[is] <- f[1]
205 vest[is] <- f[2]
206 } # end looks
207 } # end param
208 } # end sim
209
210 if(progress) cat('\n')
211
212 data.frame(sim=sim, parameter=parm, look=look,
213 est=est, vest=vest)
214 }
6767 xlab(xlab) + ylab(ylab) +
6868 guides(size = guide_legend(title='Frequency'))
6969 else
70 ggplot(k, aes(x=x, y=y, alpha=Freq ^ 0.25, label=Freq,
70 ggplot(k, aes(x=x, y=y, label=Freq, #alpha=Freq ^ 0.25,
7171 color=Freq ^ 0.25)) +
7272 geom_point(...) +
7373 scale_color_gradientn(colors=fcolors) +
8989 xlab(xlab) + ylab(ylab) +
9090 guides(size = guide_legend(title='Frequency'))
9191 else
92 ggplot(k, aes(x=x, y=y, alpha=fg, label=Freq,
92 ggplot(k, aes(x=x, y=y, label=Freq, # alpha=fg,
9393 color=if(few) k$Freq else as.integer(fg))) +
9494 geom_point(...) +
9595 scale_color_gradientn(colors=fcolors,
0 ##' Pair-up and Compute Differences
1 ##'
2 ##' This function sets up for plotting half-width confidence intervals for differences, sorting by descending order of differences within major categories, especially for dot charts as produced by [dotchartpl()]. Given a numeric vector `x` and a grouping (superpositioning) vector `group` with exactly two levels, computes differences in possibly transformed `x` between levels of `group` for the two observations that are equal on `major` and `minor`. If `lower` and `upper` are specified, using `conf.int` and approximate normality on the transformed scale to backsolve for the standard errors of estimates, and uses approximate normality to get confidence intervals on differences by taking the square root of the sum of squares of the two standard errors. Coordinates for plotting half-width confidence intervals are also computed. These intervals may be plotted on the same scale as `x`, having the property that they overlap the two `x` values if and only if there is no "significant" difference at the `conf.int` level.
3 ##' @title pairUpDiff
4 ##' @param x a numeric vector
5 ##' @param major an optional factor or character vector
6 ##' @param minor an optional factor or character vector
7 ##' @param group a required factor or character vector with two levels
8 ##' @param refgroup a character string specifying which level of `group` is to be subtracted
9 ##' @param lower an optional numeric vector giving the lower `conf.int` confidence limit for `x`
10 ##' @param upper similar to `lower` but for the upper limit
11 ##' @param minkeep the minimum value of `x` required to keep the observation. An observation is kept if either `group` has `x` exceeding or equalling `minkeep`. Default is to keep all observations.
12 ##' @param sortdiff set to `FALSE` to avoid sorting observations by descending between-`group` differences
13 ##' @param conf.int confidence level; must have been the value used to compute `lower` and `upper` if they are provided
14 ##' @return a list of two objects both sorted by descending values of differences in `x`. The `X` object is a data frame that contains the original variables sorted by descending differences across `group` and in addition a variable `subscripts` denoting the subscripts of original observations with possible re-sorting and dropping depending on `sortdiff` and `minkeep`. The `D` data frame contains sorted differences (`diff`), `major`, `minor`, `sd` of difference, `lower` and `upper` confidence limits for the difference, `mid`, the midpoint of the two `x` values involved in the difference, `lowermid`, the midpoint minus 1/2 the width of the confidence interval, and `uppermid`, the midpoint plus 1/2 the width of the confidence interval. Another element returned is `dropped` which is a vector of `major` / `minor` combinations dropped due to `minkeep`.
15 ##' @author Frank Harrell
16 ##' @export
17 ##' @md
18 ##' @examples
19 ##' x <- c(1, 4, 7, 2, 5, 3, 6)
20 ##' pairUpDiff(x, c(rep('A', 4), rep('B', 3)),
21 ##' c('u','u','v','v','z','z','q'),
22 ##' c('a','b','a','b','a','b','a'), 'a', x-.1, x+.1)
23
24 pairUpDiff <- function(x, major=NULL, minor=NULL, group, refgroup,
25 lower=NULL, upper=NULL, minkeep=NULL,
26 sortdiff=TRUE, conf.int=0.95) {
27
28 n <- length(x)
29 major <- if(! length(major)) rep(' ', n) else as.character(major)
30 minor <- if(! length(minor)) rep(' ', n) else as.character(minor)
31 ## Note: R will not let you use z[cbind(...)] if one of the dimensions
32 ## is of length 1 and equal to ''; needed to use ' '
33 group <- as.character(group)
34 glev <- unique(group)
35 if(length(glev) != 2) stop('group must have two distinct values')
36 if(refgroup %nin% glev) stop('refgroup must be one of the group values')
37 altgroup <- setdiff(glev, refgroup)
38
39 mm <- c(major, minor)
40 sep <- if(! any(grepl(':', mm))) ':'
41 else
42 if(! any(grepl('|', mm))) '|'
43 else
44 if(! any(grepl(';', mm))) ';'
45 else
46 if(! any(grepl('!', mm))) '!'
47 else
48 stop('major or minor contain all delimiters :|;!')
49
50 m <- paste0(major, sep, minor)
51 u <- unique(m)
52 lu <- length(u)
53
54 lowup <- length(lower) * length(upper) > 0
55 zcrit <- qnorm((1 + conf.int) / 2)
56
57 ## See if any observations should be dropped
58 dropped <- NULL
59 if(length(minkeep)) {
60 xa <- xb <- structure(rep(NA, lu), names=u)
61 j <- group == refgroup
62 xa[m[j]] <- x[j]
63 j <- group == altgroup
64 xb[m[j]] <- x[j]
65 j <- ! is.na(xa + xb) & (xa < minkeep & xb < minkeep)
66 if(any(j)) {
67 dropped <- names(xa)[j]
68 u <- setdiff(u, dropped)
69 lu <- length(u)
70 }
71 }
72
73 xa <- xb <- sda <- sdb <- diff <- mid <- dsd <- dlower <-
74 dupper <- dlowermid <- duppermid <- structure(rep(NA, lu), names=u)
75 j <- (group == refgroup) & (m %nin% dropped)
76 w <- m[j]
77 xa[w] <- x[j]
78 if(lowup) sda[w] <- 0.5 * (upper[j] - lower[j]) / zcrit
79
80 j <- (group == altgroup) & (m %nin% dropped)
81 w <- m[j]
82 xb[w] <- x[j]
83 if(lowup) sdb[w] <- 0.5 * (upper[j] - lower[j]) / zcrit
84
85 diff[u] <- xb[u] - xa[u]
86 if(lowup) {
87 dsd[u] <- sqrt(sda[u] ^ 2 + sdb[u] ^ 2)
88 dlower[u] <- diff[u] - zcrit * dsd[u]
89 dupper[u] <- diff[u] + zcrit * dsd[u]
90 }
91 mid[u] <- (xa[u] + xb[u]) / 2.
92 if(lowup) {
93 dlowermid[u] <- mid[u] - 0.5 * zcrit * dsd[u]
94 duppermid[u] <- mid[u] + 0.5 * zcrit * dsd[u]
95 }
96
97 k <- strsplit(u, sep)
98 ma <- sapply(k, function(x) x[[1]])
99 mi <- sapply(k, function(x) x[[2]])
100
101 ww <- list(x, major, minor, group, lower, upper, subscripts=1:length(x))
102
103 Z <- if(lowup)
104 data.frame(x, major, minor, group, lower, upper,
105 subscripts=1 : length(x))
106 else
107 data.frame(x, major, minor, group, subscripts=1 : length(x))
108
109 if(length(dropped)) Z <- Z[m %nin% dropped, ]
110 if(sortdiff) {
111 m <- paste0(Z$major, sep, Z$minor)
112 ## diff[m] is a table lookup; difference will be same for both groups
113 j <- order(Z$major, ifelse(is.na(diff[m]), -Inf, - diff[m]))
114 Z <- Z[j, ]
115 }
116
117 ## Variables referenced below have already had observations dropped
118 ## due to minkeep
119 D <- data.frame(diff=diff[u], major=ma, minor=mi, sd=dsd[u],
120 lower=dlower[u], upper=dupper[u],
121 mid=mid[u], lowermid=dlowermid[u], uppermid=duppermid[u])
122 if(sortdiff) {
123 j <- order(D$major, ifelse(is.na(D$diff), -Inf, -D$diff))
124 D <- D[j, ]
125 }
126 list(X=Z, D=D, dropped=dropped)
127 }
6161 if(length(x) && (length(x) != length(p)))
6262 stop('p and x must have same length')
6363 if(length(x) && any(diff(x) <= 0)) stop('x is not sorted or has duplicates')
64
6465 if(abs(sum(p) - 1) > .00001)
6566 stop('probabilities do not sum to 1')
6667
152153 x <- d[[v[2]]]
153154 xl <- label(x, default=v[2])
154155 s <- sn <- NULL
156 sl <- NULL
155157 if(length(v) > 2) {
156158 if(length(odds.ratio))
157 stop('odds ratio may not be specified when a stratification variable is include')
159 stop('odds ratio may not be specified when a stratification variable is included')
158160 s <- d[[v[3]]]
159161 sl <- label(s, default=v[3])
160162 s <- as.factor(s)
161163 sn <- 's'
162164 }
163165 names(d) <- c('y', 'x', sn)
166 ## ggplot2 bar chart puts first category at the top
167 ## Let's put it at the bottom
168 d$y <- factor(d$y, levels=rev(levels(as.factor(d$y))))
169
170
164171 # For each x compute the vector of proportions of y categories
165172 # Assume levels are in order
173 # Put numerators and denominators into a character string so that
174 # data.table [ operator can operate on one variable
175 # The delimiter is a single space
166176 g <- function(y) {
167177 tab <- table(y)
168 tab / sum(tab)
169 }
178 structure(paste(tab, rep(sum(tab), length(tab))), names=names(tab))
179 }
180 atxt <- function(d, strat=NULL, or=FALSE) {
181 if(! or) {
182 z <- d$prop
183 num <- as.numeric(sub(' .*', '', z)) # get first number
184 denom <- as.numeric(sub('.* ', '', z)) # get second number
185 d$prop <- num / denom
186 }
187 d$txt <- paste0(sl, if(length(strat)) ': ', as.character(strat),
188 if(length(strat)) '<br>',
189 xl, ': ', as.character(d$x), '<br>',
190 yl, ': ', as.character(d$y), '<br>',
191 'Proportion: ', round(d$prop, 3),
192 if(! or) '\u2003',
193 if(! or) markupSpecs$html$frac(num, denom, size=90))
194 d
195 }
196
170197 d <- data.table(d)
171 if(! length(s)) {
198
199 if(! length(s)) { # stratification variable not present
172200 p <- d[, as.list(g(y)), by=x]
173201 pm <- melt(p, id=1, variable.name='y', value.name='prop')
202 pm <- atxt(pm)
203
174204 plegend <- guides(fill=guide_legend(title=yl))
175205 if(! length(odds.ratio)) {
176 gg <- ggplot(pm, aes(x=as.factor(x), y=prop, fill=factor(y))) +
206 gg <- ggplot(pm, aes(x=as.factor(x), y=prop, fill=factor(y), text=txt)) +
177207 geom_col() + plegend + xlab(xl) + ylab('Proportion')
178208 return(gg)
179209 }
180 } else {
181 p <- d[, as.list(g(y)), by=.(x,s)]
210 } else { # stratification variable present; odds ratio must be absent
211 p <- d[, as.list(g(y)), by=.(x, s)]
182212 pm <- melt(p, id=c('x', 's'), variable.name='y', value.name='prop')
213 pm <- atxt(pm, pm$s)
183214 plegend <- guides(fill=guide_legend(title=yl))
184 gg <- ggplot(pm, aes(x=as.factor(x), y=prop, fill=factor(y))) +
215 gg <- ggplot(pm, aes(x=as.factor(x), y=prop, fill=factor(y), text=txt)) +
185216 facet_wrap(~ s, ncol=ncol, nrow=nrow) +
186217 geom_col() + plegend + xlab(xl) + ylab('Proportion')
187218 return(gg)
188219 }
189220
221 ## odds.ratio is present
222
190223 if(! length(ref)) ref <- p$x[1]
191 propfirstx <- as.matrix(p[x == ref, -1])
224 pfx <- as.matrix(p[x == ref, -1])
225 propfirstx <- as.numeric(sub(' .*', '', pfx)) /
226 as.numeric(sub('.* ', '', pfx))
227
192228 g <- function(x) {
193229 w <- pomodm(p=propfirstx, odds.ratio=odds.ratio(x))
194230 names(w) <- levels(y)
196232 }
197233 pa <- d[, as.list(g(x)), by=x]
198234 pma <- melt(pa, id=1, variable.name='y', value.name='prop')
235 pma <- atxt(pma, or=TRUE)
199236 w <- rbind(cbind(type=1, pm),
200237 cbind(type=2, pma))
201238 w$type <- factor(w$type, 1 : 2,
202239 c('Observed', 'Asssuming Proportional Odds'))
203 ggplot(w, aes(x=as.factor(x), y=prop, fill=factor(y))) + geom_col() +
240 ggplot(w, aes(x=as.factor(x), y=prop, fill=factor(y), text=txt)) +
241 geom_col() +
204242 facet_wrap(~ type, nrow=2) +
205243 plegend + xlab(xl) + ylab('Proportion')
206244 }
307345 guides(fill=guide_legend(override.aes=list(shape=NA), order=1)) +
308346 coord_flip() + labs(y=xlab, x='')
309347 }
348
349 utils::globalVariables('txt')
2424 nam <- dimnames(x)[[2]]
2525 dimnames(h) <- list(nam, nam)
2626 dimnames(npair) <- list(nam, nam)
27 P <- matrix(2 * (1 - pt(abs(h) * sqrt(npair - 2) / sqrt(1 - h * h),
27 P <- matrix(2 * (1 - pt(abs(h) * sqrt(npair - 2) / sqrt(pmax(1 - h * h, 0.)),
2828 npair - 2)), ncol=p)
2929 P[abs(h) == 1] <- 0
3030 diag(P) <- NA
368368
369369
370370 histSpike <-
371 function(x, side=1, nint=100, frac=.05, minf=NULL, mult.width=1,
371 function(x, side=1, nint=100, bins=NULL, frac=.05, minf=NULL, mult.width=1,
372372 type=c('proportion','count','density'),
373373 xlim=range(x),
374 ylim=c(0,max(f)), xlab=deparse(substitute(x)),
374 ylim=c(0, max(f)), xlab=deparse(substitute(x)),
375375 ylab=switch(type, proportion='Proportion',
376376 count ='Frequency',
377377 density ='Density'),
378 y=NULL, curve=NULL, add=FALSE,
378 y=NULL, curve=NULL, add=FALSE, minimal=FALSE,
379379 bottom.align=type == 'density',
380380 col=par('col'), lwd=par('lwd'), grid=FALSE, ...)
381381 {
382382 type <- match.arg(type)
383 if(! add && side != 1)
384 stop('side must be 1 if add=FALSE')
385
386 if(add && type == 'count')
387 warning('type="count" is ignored if add=TRUE')
383
384 if(! minimal) {
385 if(! add && side != 1)
386 stop('side must be 1 if add=FALSE')
387
388 if(add && type == 'count')
389 warning('type="count" is ignored if add=TRUE')
390 }
388391
389392 if(length(y) > 1) {
390393 if(length(y) != length(x))
393396 if(length(curve))
394397 warning('curve ignored when y specified')
395398
396 i <- !is.na(x + y)
399 i <- ! is.na(x + y)
397400 curve <- list(x=x[i], y=y[i])
398401 }
399402
409412 f <- table(x)
410413 x <- as.numeric(names(f))
411414 } else {
412 ncut <- nint+1
413 bins <- seq(xlim[1], xlim[2], length = ncut)
414 delta <- (bins[2]-bins[1]) / 2
415 ncut <- nint + 1
416 if(! length(bins)) bins <- seq(xlim[1], xlim[2], length = ncut)
417 delta <- (bins[2] - bins[1]) / 2
415418 f <- table(cut(x, c(bins[1] - delta, bins)))
416419
417420 x <- bins
429432 f <- den$y
430433 }
431434
432 if(!add) {
435 if(! minimal && ! add) {
433436 if(grid)
434437 stop('add=FALSE not implemented for lattice')
435438
436439 plot(0, 0, xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab, type='n')
437440 }
438
441
439442 if(side == 1 || side == 3) {
440443 l <- 1:2;
441444 ax <- 1
447450 f <- f / max(f)
448451 if(length(minf)) f <- pmax(f, minf)
449452
453 if(minimal) {
454 plot.new() # sets (0,0) to (1,1)
455 del <- diff(xlim)
456 usr <- c(xlim[1] - 0.04 * del, xlim[2] + 0.04 * del, 0, 1)
457 par(mar=c(3.5, 0.5, 0.5, 0.5), mgp=c(2, 0.5, 0), usr=usr)
458 axis(1, line=0.25)
459 title(xlab=xlab)
460 segments(x, rep(0, length(x)), x, f, lwd=lwd, col=col)
461 return(invisible(xlim))
462 }
463
450464 pr <- parGrid(grid)
451465 usr <- pr$usr;
452466 pin <- pr$pin;
4747 }
4848 gg <- function(x) if(is.character(x) || is.factor(x))
4949 'categorical' else 'numeric'
50 xlabels <- sapply(X, label)
50 ## For some reason sapply is doubling names e.g. sbp.sbp
51 sapply2 <- function(data, ...) {
52 s <- sapply(data, ...)
53 names(s) <- names(data)
54 s
55 }
56 xlabels <- sapply2(X, label)
5157 xlabels <- ifelse(xlabels == '', names(xlabels), xlabels)
52 ylabels <- sapply(Y, label)
58 ylabels <- sapply2(Y, label)
5359 ylabels <- ifelse(ylabels == '', names(ylabels), ylabels)
5460
5561 for(n in names(w)) # takes care of R change stringsAsFactors=FALSE
5763
5864 structure(w, class=c('summaryS', 'data.frame'),
5965 formula=formula, fun=fun,
60 xnames=names(X), xlabels=xlabels, xunits=sapply(X, units),
66 xnames=names(X), xlabels=xlabels, xunits=sapply2(X, units),
6167 xtype=sapply(X, gg),
62 ynames=namY, ylabels=ylabels, yunits=sapply(Y, units),
68 ynames=namY, ylabels=ylabels, yunits=sapply2(Y, units),
6369 ylim=ylim, funlabel=funlabel)
6470 }
6571
411417
412418 yvarlev <- NULL
413419 for(v in levels(X$yvar)) {
414 un <- yunits[v]
420 un <- yunits[v]
415421 l <- if(ylabels[v] == v && un == '') v else
416422 labelPlotmath(ylabels[v], un, html=TRUE)
417423 yvarlev <- c(yvarlev, l)
3030 These have the property of intersecting the two proportions if and
3131 only if there is no significant difference at the \code{1 - conf.int}
3232 level.
33
34 Specify \code{fun=exp} and \code{ifun=log} if estimates and confidence
35 limits are on the log scale. Make sure that zeros were prevented in
36 the original calculations. For exponential hazard rates this can be
37 accomplished by replacing event counts of 0 with 0.5.
3338 }
3439 \usage{
3540 dotchartpl(x, major=NULL, minor=NULL, group=NULL, mult=NULL,
3641 big=NULL, htext=NULL, num=NULL, denom=NULL,
42 numlabel='', denomlabel='',
43 fun=function(x) x, ifun=function(x) x, op='-',
3744 lower=NULL, upper=NULL,
3845 refgroup=NULL, sortdiff=TRUE, conf.int=0.95,
3946 minkeep=NULL, xlim=NULL, xlab='Proportion',
4047 tracename=NULL, limitstracename='Limits',
4148 nonbigtracename='Stratified Estimates',
42 width=800, height=NULL,
49 dec=3, width=800, height=NULL,
4350 col=colorspace::rainbow_hcl)
4451 }
4552 \arguments{
6067 \code{num} is given, \code{x} is automatically added to hover text,
6168 rounded to 3 digits after the decimal point.}
6269 \item{denom}{like \code{num} but for denominators}
70 \item{numlabel}{character string to put to the right of the numerator
71 in hover text}
72 \item{denomlabel}{character string to put to the right of the
73 denominator in hover text}
74 \item{fun}{a transformation to make when printing estimates. For
75 example, one may specify \code{fun=exp} to anti-log estimates and
76 confidence limites that were computed on a log basis}
77 \item{ifun}{inverse transformation of \code{fun}}
78 \item{op}{set to for example \code{'/'} when \code{fun=exp} and
79 effects are computed as ratios instead of differences. This is used
80 in hover text.}
6381 \item{lower}{lower limits for optional error bars}
6482 \item{upper}{upper limits for optional error bars}
6583 \item{refgroup}{if \code{group} is specified and there are exactly two
87105 \item{col}{a function or vector of colors to assign to \code{group}.
88106 If a function it will be evaluated with an argument equal to the
89107 number of distinct groups.}
108 \item{dec}{number of places to the right of the decimal place for
109 formatting numeric quantities in hover text}
90110 \item{width}{width of plot in pixels}
91111 \item{height}{height of plot in pixels; computed from number of strata
92112 by default}
121141 # state/region/group
122142 d <- subset(d, city == 0)
123143 with(d,
124 dotchartpl(x, major, minor, group, num=num, denom=denom, refgroup='Male')
144 dotchartpl(x, major, minor, group, num=num, denom=denom,
145 lower=lower, upper=upper, refgroup='Male')
125146 )
126147
127148 n <- 500
154175 big=region == 'All Regions', num=freq, denom=denom)
155176 )
156177 # Note these plots can be created by plot.summaryP when options(grType='plotly')
178
179 # Plot hazard rates and ratios with confidence limits, on log scale
180 d <- data.frame(tx=c('a', 'a', 'b', 'b'),
181 event=c('MI', 'stroke', 'MI', 'stroke'),
182 count=c(10, 5, 5, 2),
183 exposure=c(1000, 1000, 900, 900))
184 # There were no zero event counts in this dataset. In general we
185 # want to handle that, hence the 0.5 below
186 d <- upData(d, hazard = pmax(0.5, count) / exposure,
187 selog = sqrt(1. / pmax(0.5, count)),
188 lower = log(hazard) - 1.96 * selog,
189 upper = log(hazard) + 1.96 * selog)
190 with(d,
191 dotchartpl(log(hazard), minor=event, group=tx, num=count, denom=exposure,
192 lower=lower, upper=upper,
193 fun=exp, ifun=log, op='/',
194 numlabel='events', denomlabel='years',
195 refgroup='a', xlab='Events Per Person-Year')
196 )
157197 }
158198 }
159199 \keyword{hplot}
0 % Generated by roxygen2: do not edit by hand
1 % Please edit documentation in R/gbayesSeqSim.r
2 \name{estSeqSim}
3 \alias{estSeqSim}
4 \title{estSeqSim}
5 \usage{
6 estSeqSim(parameter, looks, gendat, fitter, nsim = 1, progress = FALSE)
7 }
8 \arguments{
9 \item{parameter}{vector of true parameter (effects; group differences) values}
10
11 \item{looks}{integer vector of observation numbers at which posterior probabilities are computed}
12
13 \item{gendat}{a function of three arguments: true parameter value (scalar), sample size for first group, sample size for second group}
14
15 \item{fitter}{a function of two arguments: 0/1 group indicator vector and the dependent variable vector}
16
17 \item{nsim}{number of simulations (default is 1)}
18
19 \item{progress}{set to \code{TRUE} to send current iteration number to the console}
20 }
21 \value{
22 a data frame with number of rows equal to the product of \code{nsim}, the length of \code{looks}, and the length of \code{parameter}.
23 }
24 \description{
25 Simulate Comparisons For Use in Sequential Clinical Trial Simulations
26 }
27 \details{
28 Simulates sequential clinical trials. Looks are done sequentially at observation numbers given in the vector \code{looks} with the earliest possible look being at observation 2. For each true effect parameter value, simulation, and at each look, runs a function to compute the estimate of the parameter of interest along with its variance. For each simulation, data are first simulated for the last look, and these data are sequentially revealed for earlier looks. The user provides a function \code{gendat} that given a true effect of \code{parameter} and the two sample sizes (for treatment groups 1 and 2) returns a list with vectors \code{y1} and \code{y2} containing simulated data. The user also provides a function \code{fitter} with arguments \code{x} (group indicator 0/1) and \code{y} (response variable) that returns a 2-vector containing the effect estimate and its variance. \code{parameter} is usually on the scale of a regression coefficient, e.g., a log odds ratio.
29 }
30 \examples{
31 # Run 100 simulations, 5 looks, 2 true parameter values
32 # Total simulation time: 2s
33 lfit <- function(x, y) {
34 f <- rms::lrm.fit(x, y)
35 k <- length(coef(f))
36 c(coef(f)[k], vcov(f)[k, k])
37 }
38 gdat <- function(beta, n1, n2) {
39 # Cell probabilities for a 7-category ordinal outcome for the control group
40 p <- c(2, 1, 2, 7, 8, 38, 42) / 100
41
42 # Compute cell probabilities for the treated group
43 p2 <- pomodm(p=p, odds.ratio=exp(beta))
44 y1 <- sample(1 : 7, n1, p, replace=TRUE)
45 y2 <- sample(1 : 7, n2, p2, replace=TRUE)
46 list(y1=y1, y2=y2)
47 }
48
49 set.seed(1)
50 est <- estSeqSim(c(0, log(0.7)), looks=c(50, 75, 95, 100, 200),
51 gendat=gdat,
52 fitter=lfit, nsim=100)
53 head(est)
54 }
55 \seealso{
56 \code{gbayesSeqSim()}
57 }
58 \author{
59 Frank Harrell
60 }
257257 \cr
258258 \email{fh@fharrell.com}
259259 }
260 \seealso{\code{\link{gbayesSeqSim}}}
260261 \references{
261262 Spiegelhalter DJ, Freedman LS, Parmar MKB (1994): Bayesian approaches to
262263 randomized trials. JRSS A 157:357--416. Results for \code{gbayes} are derived from
0 % Generated by roxygen2: do not edit by hand
1 % Please edit documentation in R/gbayesSeqSim.r
2 \name{gbayesSeqSim}
3 \alias{gbayesSeqSim}
4 \title{gbayesSeqSim}
5 \usage{
6 gbayesSeqSim(est, asserts)
7 }
8 \arguments{
9 \item{est}{data frame created by \code{estSeqSim()}}
10
11 \item{asserts}{list of lists. The first element of each list is the user-specified name for each assertion/prior combination, e.g., \code{"efficacy"}. The other elements are, in order, a character string equal to "<", ">", or "in", a parameter value \code{cutoff} (for "<" and ">") or a 2-vector specifying an interval for "in", and either a prior distribution mean and standard deviation named \code{mu} and \code{sigma} respectively, or a parameter value (\code{"cutprior"}) and tail area \code{"tailprob"}. If the latter is used, \code{mu} is assumed to be zero and \code{sigma} is solved for such that P(parameter > 'cutprior') = P(parameter < - 'cutprior') = \code{tailprob}.}
12 }
13 \value{
14 a data frame with number of rows equal to that of \code{est} with a number of new columns equal to the number of assertions added. The new columns are named \code{p1}, \code{p2}, \code{p3}, ... (posterior probabilities), \code{mean1}, \code{mean2}, ... (posterior means), and \code{sd1}, \code{sd2}, ... (posterior standard deviations). The returned data frame also has an attribute \code{asserts} added which is the original \code{asserts} augmented with any derived \code{mu} and \code{sigma} and converted to a data frame, and another attribute \code{alabels} which is a named vector used to map \code{p1}, \code{p2}, ... to the user-provided labels in \code{asserts}.
15 }
16 \description{
17 Simulate Bayesian Sequential Treatment Comparisons Using a Gaussian Model
18 }
19 \details{
20 Simulate a sequential trial under a Gaussian model for parameter estimates, and Gaussian priors using simulated estimates and variances returned by \code{estSeqSim}. For each row of the data frame \code{est} and for each prior/assertion combination, computes the posterior probability of the assertion.
21 }
22 \examples{
23 \dontrun{
24 # Simulate Bayesian operating characteristics for an unadjusted
25 # proportional odds comparison (Wilcoxon test)
26 # For 100 simulations, 5 looks, 2 true parameter values, and
27 # 2 assertion/prior combinations, compute the posterior probability
28 # Use a low-level logistic regression call to speed up simuluations
29 # Use data.table to compute various summary measures
30 # Total simulation time: 2s
31 lfit <- function(x, y) {
32 f <- rms::lrm.fit(x, y)
33 k <- length(coef(f))
34 c(coef(f)[k], vcov(f)[k, k])
35 }
36 gdat <- function(beta, n1, n2) {
37 # Cell probabilities for a 7-category ordinal outcome for the control group
38 p <- c(2, 1, 2, 7, 8, 38, 42) / 100
39
40 # Compute cell probabilities for the treated group
41 p2 <- pomodm(p=p, odds.ratio=exp(beta))
42 y1 <- sample(1 : 7, n1, p, replace=TRUE)
43 y2 <- sample(1 : 7, n2, p2, replace=TRUE)
44 list(y1=y1, y2=y2)
45 }
46
47 # Assertion 1: log(OR) < 0 under prior with prior mean 0.1 and sigma 1 on log OR scale
48 # Assertion 2: OR between 0.9 and 1/0.9 with prior mean 0 and sigma computed so that
49 # P(OR > 2) = 0.05
50 asserts <- list(list('Efficacy', '<', 0, mu=0.1, sigma=1),
51 list('Similarity', 'in', log(c(0.9, 1/0.9)),
52 cutprior=log(2), tailprob=0.05))
53
54 set.seed(1)
55 est <- estSeqSim(c(0, log(0.7)), looks=c(50, 75, 95, 100, 200),
56 gendat=gdat,
57 fitter=lfit, nsim=100)
58 z <- gbayesSeqSim(est, asserts)
59 head(z)
60 attr(z, 'asserts')
61
62 # Compute the proportion of simulations that hit targets (different target posterior
63 # probabilities for efficacy vs. similarity)
64
65 # For the efficacy assessment compute the first look at which the target
66 # was hit (set to infinity if never hit)
67 require(data.table)
68 z <- data.table(z)
69 u <- z[, .(first=min(p1 > 0.95])), by=.(parameter, sim)]
70 # Compute the proportion of simulations that ever hit the target and
71 # that hit it by the 100th subject
72 u[, .(ever=mean(first < Inf)), by=.(parameter)]
73 u[, .(by75=mean(first <= 100)), by=.(parameter)]
74 }
75 }
76 \seealso{
77 \code{gbayes()}, \code{estSeqSim()}
78 }
79 \author{
80 Frank Harrell
81 }
4242 2-vector, the first element corresponds to \code{x} and the second to
4343 \code{y}.}
4444 \item{g}{number of quantile groups to make for frequency counts. Use
45 \code{g=0} to use frequencies continuously for color and alpha
45 \code{g=0} to use frequencies continuously for color
4646 coding. This is recommended only when using \code{plotly}.}
4747 \item{cuts}{instead of using \code{g}, specify \code{cuts} to provide
4848 the vector of cuts for categorizing frequencies for assignment to colors}
5555 \item{xlab,ylab}{axis labels. If not specified and variable has a
5656 \code{label}, thatu label will be used.}
5757 \item{fcolors}{\code{colors} argument to pass to
58 \code{scale_color_gradientn} to color code frequencies}
58 \code{scale_color_gradientn} to color code frequencies. Use
59 \code{fcolors=gray.colors(10, 0.75, 0)} to show gray
60 scale, for example. Another good choice is
61 \code{fcolors=hcl.colors(10, 'Blue-Red')}.}
5962 \item{nsize}{set to \code{TRUE} to not vary color or transparency but
6063 instead to size the symbols in relation to the number of points. Best
6164 with both \code{x} and \code{y} are discrete. \code{ggplot2}
0 % Generated by roxygen2: do not edit by hand
1 % Please edit documentation in R/pairUpDiff.r
2 \name{pairUpDiff}
3 \alias{pairUpDiff}
4 \title{pairUpDiff}
5 \usage{
6 pairUpDiff(
7 x,
8 major = NULL,
9 minor = NULL,
10 group,
11 refgroup,
12 lower = NULL,
13 upper = NULL,
14 minkeep = NULL,
15 sortdiff = TRUE,
16 conf.int = 0.95
17 )
18 }
19 \arguments{
20 \item{x}{a numeric vector}
21
22 \item{major}{an optional factor or character vector}
23
24 \item{minor}{an optional factor or character vector}
25
26 \item{group}{a required factor or character vector with two levels}
27
28 \item{refgroup}{a character string specifying which level of \code{group} is to be subtracted}
29
30 \item{lower}{an optional numeric vector giving the lower \code{conf.int} confidence limit for \code{x}}
31
32 \item{upper}{similar to \code{lower} but for the upper limit}
33
34 \item{minkeep}{the minimum value of \code{x} required to keep the observation. An observation is kept if either \code{group} has \code{x} exceeding or equalling \code{minkeep}. Default is to keep all observations.}
35
36 \item{sortdiff}{set to \code{FALSE} to avoid sorting observations by descending between-\code{group} differences}
37
38 \item{conf.int}{confidence level; must have been the value used to compute \code{lower} and \code{upper} if they are provided}
39 }
40 \value{
41 a list of two objects both sorted by descending values of differences in \code{x}. The \code{X} object is a data frame that contains the original variables sorted by descending differences across \code{group} and in addition a variable \code{subscripts} denoting the subscripts of original observations with possible re-sorting and dropping depending on \code{sortdiff} and \code{minkeep}. The \code{D} data frame contains sorted differences (\code{diff}), \code{major}, \code{minor}, \code{sd} of difference, \code{lower} and \code{upper} confidence limits for the difference, \code{mid}, the midpoint of the two \code{x} values involved in the difference, \code{lowermid}, the midpoint minus 1/2 the width of the confidence interval, and \code{uppermid}, the midpoint plus 1/2 the width of the confidence interval. Another element returned is \code{dropped} which is a vector of \code{major} / \code{minor} combinations dropped due to \code{minkeep}.
42 }
43 \description{
44 Pair-up and Compute Differences
45 }
46 \details{
47 This function sets up for plotting half-width confidence intervals for differences, sorting by descending order of differences within major categories, especially for dot charts as produced by \code{\link[=dotchartpl]{dotchartpl()}}. Given a numeric vector \code{x} and a grouping (superpositioning) vector \code{group} with exactly two levels, computes differences in possibly transformed \code{x} between levels of \code{group} for the two observations that are equal on \code{major} and \code{minor}. If \code{lower} and \code{upper} are specified, using \code{conf.int} and approximate normality on the transformed scale to backsolve for the standard errors of estimates, and uses approximate normality to get confidence intervals on differences by taking the square root of the sum of squares of the two standard errors. Coordinates for plotting half-width confidence intervals are also computed. These intervals may be plotted on the same scale as \code{x}, having the property that they overlap the two \code{x} values if and only if there is no "significant" difference at the \code{conf.int} level.
48 }
49 \examples{
50 x <- c(1, 4, 7, 2, 5, 3, 6)
51 pairUpDiff(x, c(rep('A', 4), rep('B', 3)),
52 c('u','u','v','v','z','z','q'),
53 c('a','b','a','b','a','b','a'), 'a', x-.1, x+.1)
54 }
55 \author{
56 Frank Harrell
57 }
2626 proportions stratified by a grouping variable (and optionally a stratification variable), with an optional
2727 additional graph showing what the proportions would be had proportional
2828 odds held and an odds ratio was applied to the proportions in a
29 reference group.
29 reference group. If the result is passed to \code{ggplotly}, customized
30 tooltip hover text will appear.
3031
3132 \code{propsTrans} uses \code{\link{ggplot2}} to plot all successive
3233 transition proportions. \code{formula} has the state variable on the
5556 \arguments{
5657 \item{p}{
5758 a vector of marginal cell probabilities which must add up to one.
58 The \code{i}th element specifies the probability that a patient will be
59 in response level \code{i}, averaged over the two treatment groups.
59 For \code{popower} and \code{posamsize}, The \code{i}th element specifies the probability that a patient will be
60 in response level \code{i}, averaged over the two treatment groups. For
61 \code{pomodm}, \code{p} is the vector of cell probabilities to be
62 translated under a given odds ratio.
6063 }
6164 \item{odds.ratio}{
6265 the odds ratio to be able to detect. It doesn't
6669 apply to the refernce group to get all other group's expected proportions
6770 were proportional odds to hold against the first group. Normally the
6871 function should return 1.0 when its \code{x} argument corresponds to the
69 \code{ref} group.}
72 \code{ref} group. For \code{pomodm} is the odds ratio to apply to
73 convert the given cell probabilities.}
7074 \item{n}{
7175 total sample size for \code{popower}. You must specify either \code{n} or
7276 \code{n1} and \code{n2}. If you specify \code{n}, \code{n1} and
198202 d <- expand.grid(time=1:3, reps=1:30)
199203 d$y <- sample(letters[1:5], nrow(d), replace=TRUE)
200204 propsPO(y ~ time, data=d, odds.ratio=function(time) c(1, 2, 4)[time])
205 # To show with plotly, save previous result as object p and then:
206 # plotly::ggplotly(p, tooltip='text')
201207
202208 # Add a stratification variable and don't consider an odds ratio
203209 d <- expand.grid(time=1:5, sex=c('female', 'male'), reps=1:30)
118118 # sc(a,b) means default to a if number of axes <= 3, b if >=50, use
119119 # linear interpolation within 3-50
120120
121 histSpike(x, side=1, nint=100, frac=.05, minf=NULL, mult.width=1,
121 histSpike(x, side=1, nint=100, bins=NULL, frac=.05, minf=NULL, mult.width=1,
122122 type=c('proportion','count','density'),
123123 xlim=range(x), ylim=c(0,max(f)), xlab=deparse(substitute(x)),
124124 ylab=switch(type,proportion='Proportion',
125125 count ='Frequency',
126126 density ='Density'),
127 y=NULL, curve=NULL, add=FALSE,
127 y=NULL, curve=NULL, add=FALSE, minimal=FALSE,
128128 bottom.align=type=='density', col=par('col'), lwd=par('lwd'),
129129 grid=FALSE, \dots)
130130
204204 the curve. For \code{histSpike}, interpolated \code{y} values are
205205 derived for binmidpoints.
206206 }
207 \item{minimal}{for \code{histSpike} set \code{minimal=TRUE} to draw a
208 minimalist spike histogram with no y-axis. This works best when
209 produce graphics images that are short, e.g., have a height of
210 two inches. \code{add} is forced to be \code{FALSE} in this case
211 so that a standalone graph is produced. Only base graphics are
212 used.}
207213 \item{bottom.align}{
208214 set to \code{TRUE} to have the bottoms of tick marks (for
209215 \code{side=1} or \code{side=3}) aligned at the y-coordinate. The
259265 plotted. For \code{histSpikeg}, if \code{x} has no more than
260266 \code{nint} unique values, all observed values are used, otherwise
261267 the data are rounded before tabulation so that there are no more
262 than \code{nint} intervals.
263 }
268 than \code{nint} intervals. For \code{histSpike}, \code{nint} is
269 ignored if \code{bins} is given.
270 }
271 \item{bins}{for \code{histSpike} specifies the actual cutpoints to use
272 for binning \code{x}. The default is to use \code{nint} in
273 conjunction with \code{xlim}.}
264274 \item{\dots}{
265275 optional arguments passed to \code{scat1d} from \code{datadensity}
266276 or to \code{histSpike} from \code{scat1d}. For \code{histSpikep}