Description Details Author(s) Examples
For the case of treatment-covariate interaction in linear models, this package provides functions to preform multiple comparisons among several treatments over a grid of (many) pre-specified covariate values. It thus may be used as an alternative to confidence bands for differences of multiple treatments. For the special case of (general) linear models with one treatment factor and one covariate, simultaneous confidence intervals for between-treatment-contrasts over a prespecified grid of covariate values can be computed. Such confidence intervals can be computed for differences of regression lines and (approximately) ratios of regression lines. For use in linear or generalized linear models, matrices of linear combinations for multiple pairwise between-treatment comparisons can be computed that can be used for plug-in into the linfct argument of glht, package multcomp.
Package: | statintcov |
Type: | Package |
Version: | 1.0 |
Date: | 2015-03-09 |
License: | What license is it under? |
For the special case of one treatment factor and interaction to one covariate, the function scitreatcov
has a simple interface allowing the specification of response, tretament factor and covariate, the treatment contrast (for differences) of interest and a vector of covariate values. Internally, a suitably parameterized linear model is fit, a corresponding contrast matrix is contructed. Returned are simultaneous confidence intervals for that grid, internally using the function confint.glht
from package multcomp
. The function sciratiotreatcov
allows to perform the same type of analysis for ratios of regression lines, relying on the approximate (!) methods available in function gsci.ratio
from package mratios
.
For the slightly more general case of (generalized) linear models, that may involve more than one covariate or linear and quadratic dependency on the covariate, the function cmiacov
is provided. For a given fitted model object, a grid of covariate values, and a matrix defining pairwise differences of interest among the levels of a treatment factor, a matrix of coefficients is computed that can be plugged into the argument linfct
of the glht
function in package multcomp
.
Frank Schaarschmidt <schaarschmidt@biostat.uni-hannover.de>
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 | ###############################################
# All pairwise comparisons in a linear model #
# with baseline as covariate #
###############################################
# Loading the data set
if(require("MASS")){
# load data from package MASS:
data(anorexia, package="MASS")
ggplot(anorexia, aes(y=Postwt, x=Prewt)) + geom_point() +
facet_wrap(~Treat)
fitan <- lm(Postwt ~ Treat*Prewt, data=anorexia)
anova(fitan)
# Comparison to control (cont, which is the second level of factor Treat)
# does any of the two treatments leads to higher postweights compared to control
# at covariate values 75,76,...,95
dscian <- scitreatcov(response="Postwt", treatment="Treat",
covariate="Prewt", data=anorexia, covset=seq(from=75, to=95, by=5),
treatcon="Dunnett", conf.level=0.95, alternative="greater", base=2)
str(dscian, max.level=1)
str(dscian$sci)
ggplot(dscian$sci, aes(x=covariate, y=Estimate, ymin=lwr, ymax=upr)) +
geom_errorbar(width=1) + geom_line() + geom_point(shape=15) +
facet_grid(.~comparison) + xlab("Preweight") +
ylab("Difference in expected postweight") + geom_hline(yintercept=0)
}
###############################################
# Comparisons to control with interaction to #
# a quadratic regression term #
###############################################
data(pc)
ggplot(pc, aes(y=yield, x=x)) + geom_point() + facet_wrap(~additive)
# Model selected by Milliken and Johnson:
fitpc<-lm(yield ~ x + I(x^2) + additive + additive:I(x^2), data=pc)
anova(fitpc)
# Since this model contains a quadratic term, it can onyl be handled via cmiacov
# 1) computing the contrast matrix, suitable for parameterization in fitpc
# the only covariate in the data is x, although it appears in several terms
# comparsisons to control
cmpc <- cmiacov(fitpc, treatment="additive",
covset=list("x" = seq(from=1, to=10, by=1)), treatcon="Dunnett")
str(cmpc, max.level=1)
str(cmpc$linfct)
# 2) pass the contrast matrix cmpc$linfct to the linfct argument of glht
scipc <- confint(glht(fitpc, linfct=cmpc$linfct))
# 3) for plotting: write the comparsions and covariate values in a data.frame
# together with the confidence interavls from glht
dscipc <- cbind(cmpc$datadiff, scipc$confint)
ggplot(dscipc, aes(x=x, y=Estimate, ymin=lwr, ymax=upr)) +
geom_errorbar(width=1) + geom_line() + geom_point(shape=15) +
facet_grid(.~comparison) + xlab("Compound x") +
ylab("Difference in expected yield") + geom_hline(yintercept=0)
######################################################
# All pairwise comparisons in a binomial generalized #
# linear model with logit link #
######################################################
# to load the data set
if(require("drc")){
data(selenium, package="drc")
se <- subset(selenium, conc!=0)
se$typef <- factor(se$type)
levels(se$typef) <- c("Selenate","Selenite","Selenomethionine","Selenocysteine")
se$l10conc <- log10(se$conc)
se$mortality <- se$dead/se$total
ggplot(se, aes(y=mortality, x=conc)) + geom_point() +
facet_wrap( ~typef) + scale_x_log10()
fitse<-glm(cbind(dead, total-dead) ~ typef*l10conc, data=se, family=binomial())
anova(fitse, test="Chisq")
# the binomial assumption, or logit-log linearity may not be adequat
# using the quasibinomial assumption should be more reliable
# 1) computing the contrast matrix, suitable for parameterization in fitse
# all pairwise comparisons ("Tukey")
# a sequence of 7 points over the observed range of the covariate
range(se$l10conc)
cmse <- cmiacov(fitse, treatment="typef",
covset=list("l10conc" = seq(from=min(se$l10conc),
to=max(se$l10conc), length.out=7)), treatcon="Tukey")
str(cmse, max.level=1)
str(cmse$linfct)
# 2) pass the contrast matrix cman$linfct to the linfct argument of glht
scise <- confint(glht(fitse, linfct=cmse$linfct))
# 3) for plotting: write the comparisons and covariate values in a data.frame
# together with the confidence intervals from glht
dscise <- cbind(cmse$datadiff, scise$confint)
ggplot(dscise, aes(x=l10conc, y=Estimate, ymin=lwr, ymax=upr)) +
geom_errorbar(width=0.1) + geom_line() + geom_point(shape=15) +
facet_wrap(~comparison) + xlab("Log10 Concentration") +
ylab("Difference (logit-scale)") + geom_hline(yintercept=0)
}
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.