R Under development (unstable) (2024-01-24 r85824 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > ## install.packages("~/R/regdevelop/pkg/relevance_1.3.tar.gz", repos=NULL, lib="~/local/R_libs") > require(relevance) Loading required package: relevance > > ## > d.permeability <- + data.frame(perm = c(0.80, 0.83, 1.89, 1.04, 1.45, 1.38, 1.91, 1.64, 0.73, 1.46, + 1.15, 0.88, 0.90, 0.74, 1.21), atterm = rep(1:0, c(10,5)) + ) > rtt <-t.test(perm~atterm, data=d.permeability, var.equal=FALSE) > rr <-twosamples(perm~atterm, data=d.permeability, var.equal=FALSE) > stopifnot(all(abs(rtt$conf.int+rr[c("ciUp","ciLow")])<1e-10)) > stopifnot(abs( rtt$p.value-rr["p.value"])<1e-10) > > rr1 <- twosamples(rep(0:1,c(5,20))) > onesample(rep(0:1,c(5,20))) One proportion, binomial inference log odds : 1.386294 ; confidence int.: [ 0.3762261, 2.6129210 ] Rle: 13.863 ; Rlp: 26.129 ; Rls: 3.762 ++ Relevance codes: -Inf 0 . 1 + 2 ++ 5 +++ Inf Relevance threshold: prop = 0.1 estimate: probability of success 0.8 > rr2 <- twosamples(rep(0:1,c(5,20)), rep(c(0:1,0:1),c(2,3,12,8)), rlv.threshold=0.1) > (rrf <- fisher.test(rep(0:1,c(5,20)), rep(c(0:1,0:1),c(2,3,12,8)))) Fisher's Exact Test for Count Data data: rep(0:1, c(5, 20)) and rep(c(0:1, 0:1), c(2, 3, 12, 8)) p-value = 0.6232 alternative hypothesis: true odds ratio is not equal to 1 95 percent confidence interval: 0.03154351 4.98858689 sample estimates: odds ratio 0.4594844 > > ## one sample > data(sleep) > dd <- subset(sleep, group==2) > rr <- onesample(60*dd$extra, rlv.threshold=60, standardize=FALSE) > onesample(I(60*extra) ~ 1, data=sleep, subset=group==2, + rlv.threshold=60, standardize=FALSE) One Sample t inference target variable: I(60 * extra) mean : 139.8 ; confidence int.: [ 53.86065, 225.73935 ] Rle: 2.33 ; Rlp: 3.762 ; Rls: 0.898 . Relevance codes: -Inf 0 . 1 + 2 ++ 5 +++ Inf Relevance threshold: = 60 estimate: mean se 139.80 37.99 > twosamples(I(60*extra) ~ group, data=sleep, rlv.threshold=60, standardize=FALSE) Two Sample t inference, equal variances assumed target variable: I(60 * extra) difference of means: 94.8 ; confidence int.: [ -12.23244, 201.83244 ] Rle: 1.58 ; Rlp: 3.364 ; Rls: -0.204 Relevance codes: -Inf 0 . 1 + 2 ++ 5 +++ Inf Relevance threshold: = 60 estimate: mean se 45.00000 33.94407 139.80000 37.99000 > > topt <- options(contrasts=c("contr.sum", "contr.poly")) > rr <- lm(I(60*extra) ~ group, data=sleep) > rte <- termeffects(rr) > plot(termeffects(rr)) > options(topt) ## restore options > ## ------------------------------------------------------------ > ff <- function(dd, ...) twosamples(perm~atterm, data=dd, ...) > ff(d.permeability) Two Sample t inference, equal variances assumed target variable: perm difference of means: 0.337 ; confidence int.: [ -0.1162891, 0.7902891 ] Rle: 4.147 ; Rlp: 9.725 ; Rls: -1.431 Relevance codes: -Inf 0 . 1 + 2 ++ 5 +++ Inf Relevance threshold: stand = 0.1 estimate: mean se 0.97600000 0.08823831 1.31300000 0.13952340 > ff <- function(x, ...) twosamples(x, ...) > dd <- subset(sleep, group==2) > ff(dd$extra) One Sample t inference mean : 2.33 ; confidence int.: [ 0.8976775, 3.7623225 ] Rle: 11.637 ; Rlp: 18.79 ; Rls: 4.483 ++ Relevance codes: -Inf 0 . 1 + 2 ++ 5 +++ Inf Relevance threshold: stand = 0.1 estimate: mean se 2.3300000 0.6331666 > twosamples(dd$extra) One Sample t inference mean : 2.33 ; confidence int.: [ 0.8976775, 3.7623225 ] Rle: 11.637 ; Rlp: 18.79 ; Rls: 4.483 ++ Relevance codes: -Inf 0 . 1 + 2 ++ 5 +++ Inf Relevance threshold: stand = 0.1 estimate: mean se 2.3300000 0.6331666 > > ## ================================================================== > > data(d.blast) > rlm <- lm(log10(tremor)~location+log10(distance)+log10(charge), data=d.blast) > rt <- termtable(rlm) > rt lm : Drop-term inference data: d.blast ; target variable: log10(tremor) coef df R2x coefRlp coefRls dropRls..sy predRle (Intercept) 2.964 1 . 15.27 13.18 . . location . 7 0.102 . . 5.49 +++ 3.65 log10(distance) -1.518 1 0.477 13.63 11.53 11.42 +++ 9.48 log10(charge) 0.636 1 0.102 9.62 7.53 7.46 +++ 5.50 Relevance codes: -Inf 0 . 1 + 2 ++ 5 +++ Inf Relevance thresholds: coef = 0.1, drop = 0.1, pred = 0.05 > inference(rlm) lm data: d.blast ; target variable: log10(tremor) $ termtable coef df R2x coefRlp coefRls dropRls..sy predRle (Intercept) 2.964 1 . 15.27 13.18 . . location . 7 0.102 . . 5.49 +++ 3.65 log10(distance) -1.518 1 0.477 13.63 11.53 11.42 +++ 9.48 log10(charge) 0.636 1 0.102 9.62 7.53 7.46 +++ 5.50 $ termeffects $location loc1 loc2 loc3 loc4 loc5 0.00000 0.15306 + 0.13169 + -0.16185 + -0.03211 loc6 loc7 loc8 0.07161 . -0.00889 0.00372 Relevance codes: -Inf 0 . 1 + 2 ++ 5 +++ Inf Relevance thresholds: coef = 0.1, drop = 0.1, pred = 0.05 Multiple R^2: 0.795; Adjusted R^2: 0.79 F-statistic: 152 on 9 and 352 d.f.; p.value: 1.82e-115 > rte <- termeffects(rlm) > print(rte, show.inference=c("classical","coefRls","coefRls.symbol"), single=TRUE) lm : Term effects $ (Intercept) 2.963865 $ location loc1 loc2 loc3 loc4 loc5 0.00000 0.15306 + 0.13169 + -0.16185 + -0.03211 loc6 loc7 loc8 0.07161 . -0.00889 0.00372 $ log10(distance) -1.52 +++ $ log10(charge) 0.636 +++ Relevance codes: -Inf 0 . 1 + 2 ++ 5 +++ Inf > plot(rte, single=TRUE) > > rr <- inference(rlm) > rlm <- lm(log10(tremor)~location+log10(distance)+log10(charge), + data=d.blast, subset=location %in% c("loc2", "loc4", "loc6")) > > rpr <- print(termeffects(rlm), print=FALSE) > attr(rpr, "head") <- sub("lm", "Linear Regression", attr(rpr, "head")) > rpr Linear Regression : Term effects $ location loc2 loc4 loc6 0.0000 -0.2908 +++ 0.0051 Relevance codes: -Inf 0 . 1 + 2 ++ 5 +++ Inf > > data(swiss) > rr <- lm(Fertility ~ . , data = swiss) > rt <- termtable(rr) > rt lm : Drop-term inference data: swiss ; target variable: Fertility coef df R2x coefRlp coefRls dropRls..sy predRle (Intercept) 66.915 1 . 12.76 6.53 . . Agriculture -0.172 1 0.562 6.89 0.66 0.50 . 1.12 Examination -0.258 1 0.728 4.68 -1.55 0.00 0.01 Education -0.871 1 0.640 10.46 4.23 3.76 ++ 4.16 Catholic 0.104 1 0.484 7.67 1.44 1.30 + 1.69 Infant.Mortality 1.077 1 0.097 7.47 1.24 1.11 + 1.53 Relevance codes: -Inf 0 . 1 + 2 ++ 5 +++ Inf Relevance thresholds: coef = 0.1, drop = 0.1, pred = 0.05 > rtp <- print(rt) lm : Drop-term inference data: swiss ; target variable: Fertility coef df R2x coefRlp coefRls dropRls..sy predRle (Intercept) 66.915 1 . 12.76 6.53 . . Agriculture -0.172 1 0.562 6.89 0.66 0.50 . 1.12 Examination -0.258 1 0.728 4.68 -1.55 0.00 0.01 Education -0.871 1 0.640 10.46 4.23 3.76 ++ 4.16 Catholic 0.104 1 0.484 7.67 1.44 1.30 + 1.69 Infant.Mortality 1.077 1 0.097 7.47 1.24 1.11 + 1.53 Relevance codes: -Inf 0 . 1 + 2 ++ 5 +++ Inf Relevance thresholds: coef = 0.1, drop = 0.1, pred = 0.05 > plot(rt) > ## ---- > ## glm > ## Dobson (1990) Page 93: Randomized Controlled Trial : > d.dobson <- data.frame(treatment=gl(3,3), outcome=gl(3,1,9), + counts=c(18,17,15,20,10,20,25,13,12)) > rglm <- glm(counts ~ outcome + treatment, data=d.dobson, family = poisson()) > summary(rglm) Call: glm(formula = counts ~ outcome + treatment, family = poisson(), data = d.dobson) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 3.045e+00 1.709e-01 17.815 <2e-16 *** outcome2 -4.543e-01 2.022e-01 -2.247 0.0246 * outcome3 -2.930e-01 1.927e-01 -1.520 0.1285 treatment2 1.217e-15 2.000e-01 0.000 1.0000 treatment3 8.438e-16 2.000e-01 0.000 1.0000 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 10.5814 on 8 degrees of freedom Residual deviance: 5.1291 on 4 degrees of freedom AIC: 56.761 Number of Fisher Scoring iterations: 4 > rt <- termtable(rglm) > ## relevance:::print.inference(rt, show="test") > print(rt, show="test") glm ; family = poisson : Drop-term inference data: d.dobson ; target variable: counts coef df ciLow ciUp R2x Sig0 p.value..sy (Intercept) 3.04 1 2.57 3.52 . . . outcome . 2 . . 0 0.95 0.065 . treatment . 2 . . 0 0.00 1.000 Significance codes for p.value: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1 > (rte <- termeffects(rglm)) glmlm ; family = poisson : Term effects $ outcome 1 2 3 0.000 -0.454 -0.293 $ treatment 1 2 3 0.00e+00 1.22e-15 8.44e-16 Relevance codes: -Inf 0 . 1 + 2 ++ 5 +++ Inf > rr <- inference(rglm) > ## an example with offsets from Venables & Ripley (2002, p.189) > data(anorexia, package = "MASS") > rglma <- glm(Postwt ~ Prewt + Treat + offset(Prewt), + family = gaussian, data = anorexia) > summary(rglma) Call: glm(formula = Postwt ~ Prewt + Treat + offset(Prewt), family = gaussian, data = anorexia) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 49.7711 13.3910 3.717 0.000410 *** Prewt -0.5655 0.1612 -3.509 0.000803 *** TreatCont -4.0971 1.8935 -2.164 0.033999 * TreatFT 4.5631 2.1333 2.139 0.036035 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for gaussian family taken to be 48.69504) Null deviance: 4525.4 on 71 degrees of freedom Residual deviance: 3311.3 on 68 degrees of freedom AIC: 489.97 Number of Fisher Scoring iterations: 2 > termtable(rglma) glm ; family = gaussian : Drop-term inference data: anorexia ; target variable: Postwt coef df R2x coefRlp coefRls dropRls..sy predRle (Intercept) 49.771 1 . 6.88 2.07 . . Prewt -0.566 1 0.017 6.63 1.82 1.68 + 1.48 Treat . 2 0.009 . . 3.79 ++ 3.36 Relevance codes: -Inf 0 . 1 + 2 ++ 5 +++ Inf Relevance thresholds: coef = 0.1, drop = 0.1, pred = 0.05 > rte <- termeffects(rglma) > > ## A Gamma example, from McCullagh & Nelder (1989, pp. 300-2) > clotting <- data.frame( + u = c(5,10,15,20,30,40,60,80,100), + lot1 = c(118,58,42,35,27,25,21,19,18), + lot2 = c(69,35,26,21,18,16,13,12,12)) > rglmc1 <- glm(lot1 ~ log(u), data = clotting, family = Gamma) > rglmc2 <- glm(lot2 ~ log(u), data = clotting, family = Gamma) > > rt <- termtable(rglmc1) > ## Aliased ("S"ingular) -> 1 NA coefficient > rglmc3 <- glm(lot2 ~ log(u) + log(u^2), data = clotting, family = Gamma) > ## does not give the coefficient of log(u) > > data(housing, package="MASS") > rpolr <- MASS::polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing) > rt <- termtable(rpolr) Re-fitting to get Hessian Warning message: In termtable(rpolr) : :termtable: drop1 did not work. I return the coefficient table > rt data: housing ; target variable: Sat coef df coefRlp coefRls InflMedium 0.566 1 1.80 0.84 InflHigh 1.289 1 2.96 2.00 TypeApartment -0.572 1 1.65 0.69 TypeAtrium -0.366 1 1.06 0.10 TypeTerrace -1.091 1 2.24 1.28 ContHigh 0.360 1 1.40 0.44 Relevance codes: -Inf 0 . 1 + 2 ++ 5 +++ Inf Warning message: In print.coeftable(x, lshow, print = FALSE) : :print.inference: column dropRls not available > > ##- data(d.surveyenvir) > ##- rpolr2 <- MASS::polr(disturbance ~ age+education+location, data=d.surveyenvir) > ##- rt <- termtable(rpolr2) > > ## ---------------------- > data(ovarian) ## , package="survival" > rsrw <- survival::survreg(survival::Surv(futime, fustat) ~ ecog.ps + rx, + data=ovarian, dist='weibull', scale=1) > termtable(rsrw) survreg ; distribution = weibull : Drop-term inference data: ovarian ; target variable: survival::Surv(futime, fustat) coef df R2x coefRlp coefRls dropRls..sy predRle (Intercept) 6.962 1 . 14.97 6.53 . . ecog.ps -0.433 1 1 5.73 -2.72 0 0 rx 0.582 1 1 6.24 -2.20 0 0 Relevance codes: -Inf 0 . 1 + 2 ++ 5 +++ Inf Relevance thresholds: coef = 0.1, drop = 0.1, pred = 0.05 > summary(rsrw) Call: survival::survreg(formula = survival::Surv(futime, fustat) ~ ecog.ps + rx, data = ovarian, dist = "weibull", scale = 1) Value Std. Error z p (Intercept) 6.962 1.322 5.27 1.4e-07 ecog.ps -0.433 0.587 -0.74 0.46 rx 0.582 0.587 0.99 0.32 Scale fixed at 1 Weibull distribution Loglik(model)= -97.2 Loglik(intercept only)= -98 Chisq= 1.67 on 2 degrees of freedom, p= 0.43 Number of Newton-Raphson Iterations: 4 n= 26 > > rsre <- survival::survreg(survival::Surv(futime, fustat) ~ ecog.ps + rx, + data=ovarian, dist="exponential") > termtable(rsre) survreg ; distribution = exponential : Drop-term inference data: ovarian ; target variable: survival::Surv(futime, fustat) coef df R2x coefRlp coefRls dropRls..sy predRle (Intercept) 6.962 1 . 14.97 6.53 . . ecog.ps -0.433 1 1 5.73 -2.72 0 0 rx 0.582 1 1 6.24 -2.20 0 0 Relevance codes: -Inf 0 . 1 + 2 ++ 5 +++ Inf Relevance thresholds: coef = 0.1, drop = 0.1, pred = 0.05 > > data(tobin, package="survival") > rsrtobin <- + survival::survreg(survival::Surv(durable, durable>0, type='left') ~ age + quant, + data=tobin, dist='gaussian') > rt <- termtable(rsrtobin) > > proc.time() user system elapsed 1.20 0.12 1.32