Simultaneous confidence intervals for Simpson indices

Description

NOTE: This is a Test-version and is not sufficiently checked for correctness so far. Simultaneous confidence intervals for differences and ratios of Simpsons indices of diversity are calculated from data sets with repeated samples of communities and designs with more than two treatments groups. The intervals are calculated based on a stratified, non-parametric ordinary bootstrap sample of Simpsonindices, and applying the Algorithm of Besag et al.(1995) on the joint empirical distribution of differences (BOOTSimpsonD) or ratios (BOOTSimpsonR) of the original distribution.

Usage

BOOTSimpsonD(X, f, type = "Dunnett", cmat = NULL, conf.level = 0.95,
 alternative=c("two.sided", "less", "greater"), madj=TRUE, ...)

BOOTSimpsonR(X, f, type = "Dunnett", cmat = NULL, conf.level = 0.95,
 alternative=c("two.sided", "less", "greater"), madj=TRUE, ...)

Arguments

Details

NOTE: This is a test version!

Value

If madj=TRUE, an object of class "SCSnp", see SCSnp for details. If madj=FALSE, an object of class "CInp", see CInp for details.

Author(s)

Frank Schaarschmidt

See Also

SCSnp, these function internally make use of CCDiff.boot, CCDiff.default, CCRatio.boot, CCRatio.default, boot and estSimpsonf.

Examples

X<-t(rmultinom(n=40,size=100,
 prob=c(0.3,0.2,0.2,0.1,0.1,0.05,0.05)))

colnames(X)<-paste("Sp",1:7,sep="")

DAT<-as.data.frame(X)

f<-as.factor(rep(c("A","B","C","D"),each=10))

SCIdunnettd<-BOOTSimpsonD(X=DAT, f=f, type = "Dunnett",
 conf.level = 0.95, alternative = "two.sided", R=500)
 # with small no of bootstraps for illustration

SCIdunnettd

Eklektor counts of Brachycera

Description

In a field trial, 4 treatments were arranged in a randomized complete block design with 8 blocks and 32 plots. Soil eklektor traps were placed in each plot, on six dates from 2005-07-12 to 2005-09-25, the number of individuals of Brachycera (Flies, Order Diptera) hatching from soil were counted. The individuals were classified to the family level. Interest was in assessing potential effects of the novel treatment (Novum) on the abundance of Brachycera, compared to a near standard (Standard) and two additional standard treatments, A and B.

Usage

data(Brachycera)

Format

A data frame with 192 observations on the following 15 variables.

Date

a POSIXt variable, the time of counting the individuals in the eklektor trap

Treatment

a factor with 4 levels A B Standard Novum, where Novum is the novel treatment of interest in safety assessment, and Standard is the nearest standard treatment which commonly accepted. A and B are two additional standard treatments.

Block

a numeric vector, specifying the eight blocks 1-8

Plot

a factor with levels A1 A2 to Standard8, indicator of the individuals plots

Agromy

a numeric vector, counts of individuals

Anthom

a numeric vector, counts of individuals

Callip

a numeric vector, counts of individuals

Chloro

a numeric vector, counts of individuals

Ephyd

a numeric vector, counts of individuals

Droso

a numeric vector, counts of individuals

Hybo

a numeric vector, counts of individuals

Musci

a numeric vector, counts of individuals

Phori

a numeric vector, counts of individuals

Sphaer

a numeric vector, counts of individuals

Total

a numeric vector, counts of individuals

Source

...

Examples

data(Brachycera)

par(mar=c(11,5,3,1))

boxplot(Total ~ Treatment*Date, data=Brachycera, las=2, 
col=c("white","white","blue","green"))

legend(x=15, y=80, legend=levels(Brachycera$Treatment),
 fill=c("white","white","blue","green"))

Contrasts of parameters simulated in BUGS

Description

Calculate linear combinations of parameters simulated in BUGS. This is a test version for internal use!

Usage

CCDiff(bugs, dat, cmat = NULL, 
type = c("Dunnett", "Tukey", "Sequen", "Williams", "Changepoint"))

Arguments

Details

Testversion, for internal use.

Value

An object of class "CCDiff", a list with elements

Compute contrasts of chains of joint empirical distributions.

Description

Compute contrasts of chains of joint empirical distributions obtained by stratified bootstrap. For internal use.

Usage

CCDiff.boot(x, cmat = NULL,
 type = c("Dunnett", "Tukey", "Sequen", "Williams",
 "Changepoint", "McDermott", "GrandMean", "Marcus"))

Arguments

Details

Testversion. For internal use.

Compute contrasts of chains of joint empirical distributions.

Description

Compute contrasts of chains of joint empirical distributions. For internal use.

Usage

CCDiff.default(x, cmat)

Arguments

Details

Denote the elements of x by x_{nk} and denote the elements of cmat by c_{mk}. Function CCDiff.default simply calculates:

\sum_{k=1}^{K}c_{mk}x_{nk}

for each m=1,...,M and n=1,...,N. The result is a N times M matrix.

Examples

# What the function does:

# a 10 times 4 matrix

X<-round(cbind(
 rnorm(10,1,1),
 rnorm(10,1,1),
 rnorm(10,1,1),
 rnorm(10,1,1)))

# and a x times 4 contrast matrix

CMAT<-rbind(
c(-1,1,0,0),
c(-1,0,1,0),
c(-1,0,0,1)
)

CCDiff.default(x=X, cmat=CMAT)

Ratio - contrasts of parameters simulated in BUGS

Description

Calculate ratios of parameters simulated in BUGS. This is a test version for internal use!

Usage

CCRatio(bugs, dat, cmat = NULL,
 type = c("Dunnett", "Tukey", "Sequen", "Williams", "Changepoint"))

Arguments

Details

A testversion.

Value

An object of class "CCRatio", a list with elements

Compute ratio contrasts of chains of joint empirical distributions.

Description

Compute contrasts of chains of joint empirical distributions. For internal use.

Usage

CCRatio.boot(x, cmat = NULL,
 type = c("Dunnett", "Tukey", "Sequen", "Williams",
 "Changepoint", "McDermott", "GrandMean", "Marcus"))

Arguments

Details

Testversion. For internal use.

Compute ratio contrasts of chains of joint empirical distributions.

Description

Compute ratio contrasts of chains of joint empirical distributions. For internal use.

Usage

CCRatio.default(x, cmat)

Arguments

Details

Denote the elements of x by x_{nk}. Denote the numetator of cmat by C with elements c_{mk} and the denominator of cmat as D with elements d_{mk}. Function CCRatio.default simply calculates

\frac{\sum_{k=1}^{K}c_{mk}x_{nk}}{\sum_{k=1}^{K}d_{mk}x_{nk}}

for each m=1,...,M and n=1,...,N. The result is a N times M matrix.

Examples

X<-round(cbind(
 rnorm(10,1,1),
 rnorm(10,1,1),
 rnorm(10,1,1),
 rnorm(10,1,1)))

# and numerator and denominator
# x times 4 contrast matrix

NMAT<-rbind(
c(1,0,0,0),
c(1,0,0,0),
c(1,0,0,0)
)


DMAT<-rbind(
c(0,1,0,0),
c(0,0,1,0),
c(0,0,0,1)
)


CCRatio.default(x=X, cmat=list(numC=NMAT, denC=DMAT) )

Wrapper to compute confidence intervals from glms

Description

Computes confidence intervals from the output of a glm, by calling to glht(multcomp).

Usage

CIGLM(x, conf.level = 0.95, method = c("Raw", "Adj", "Bonf"))

Arguments

Details

This is just a wrapper to confint.glht of package multcomp.

Note that except for the simple general linear model with assumption of Gaussian response, the resulting intervals are exact intervals. In other cases, the methods are only asymptotically correct, hence might give misleading results for small sample sizes!

Value

An object of class "confint.glht"

See Also

confint.glht in package multcomp for the function that is used internally, UnlogCI for a simple function to bring confidence intervals back to the original scales when there is a log or logit link, with appropriate naming.

Examples

data(Diptera)
library(multcomp)

modelfit <- glm(Ges ~ Treatment, data=Diptera, family=quasipoisson)
comps <- glht(modelfit, mcp(Treatment="Tukey"))
CIs<-CIGLM(comps, method="Raw")
CIs

CIsAdj<-CIGLM(comps, method="Adj")
CIsAdj

CIsBonf<-CIGLM(comps, method="Bonf")
CIsBonf

library(gamlss)

modelfit2 <- gamlss(Ges ~ Treatment, data=Diptera, family=NBI)
comps2 <- glht(modelfit2, mcp(Treatment="Tukey"))
CIs2<-CIGLM(comps2, method="Raw")
CIs2

CIsAdj2<-CIGLM(comps2, method="Adj")
CIsAdj2

CIsBonf2<-CIGLM(comps2, method="Bonf")
CIsBonf2

Construct local confidence intervals from joint empirical distribution.

Description

Construct local confidence intervals for each parameter from the empirical joint distribution of a parameter vector of length P.

Usage

## Default S3 method:
CInp(x, conf.level = 0.95,
 alternative = "two.sided", ...)

## S3 method for class 'CCRatio'
CInp(x, ...)

## S3 method for class 'CCDiff'
CInp(x, ...)

## S3 method for class 'bugs'
CInp(x, conf.level = 0.95,
 alternative = "two.sided", whichp = NULL, ...)

Arguments

Details

Construct simple confidence intervals based on order statistics applied to the marginal empirical distributions in x.

Value

An object of class "CInp", a list with elements

See Also

The function internally used is quantile with its default settings. See SCSnp for simultaneous sets.

Examples

# Assume a 100 times 4 matrix of 4 mutually independent
# normal variables:

X<-cbind(rnorm(100), rnorm(100), rnorm(100), rnorm(100))

lcits<-CInp(x=X, conf.level=0.95, alternative="two.sided")
lcits

ci1<-lcits$conf.int[1,]
length( which(X[,1]>=ci1[1] & X[,1]<=ci1[2] ) )


ci2<-lcits$conf.int[2,]
length( which(X[,2]>=ci2[1] & X[,2]<=ci2[2] ) )

Catches of Planthoppers and Leafhoppers

Description

Data of a field trial concerning the impact of a genetically modified variety on the abundance of Planthoppers and Leafhoppers. The trial was designed as a randomized complete block design with 8 blocks (Row). In each block, three treatments were randomized: a conventional variety treated with insecticides (Insecticide), a genetically modified variety (GM), and the near-isogenic line (Iso) the to genetically modified line.

Usage

data(Cica1)

Format

A data frame with 24 observations on the following 6 variables.

Field

a factor with levels 1 2, separating the two major sites of the trial. On field 1, the blocks 1-5 were situated, on field 2, blocks 6-8 were situated.

Row

a factor with 8 levels, specifying the blocks:R1 R2 R3 R4 R5 R6 R7 R8

Year

a numeric vector, for year 1 of the trial

Treatment

a factor with 3 levels, specifying the genetically modified variety GM, the conventional variety treated with insecticides Insecticide, and the variety that was near-isogenic to the GM variety Iso

Au_Bonitur

Counts of Auchenorryhncha by visual assessment

Zs_sweep_netting

Counts of the major species Zyginidia scutellaris, catched by sweep nets

Source

...

Examples

data(Cica1)

layout(matrix(1:2,ncol=1))

ylim<-range(Cica1[,c("Au_Bonitur","Zs_sweep_netting")])

boxplot(Au_Bonitur ~ Treatment, data=Cica1,
 main= "Aucherrhyncha, visual assessment", ylim=ylim)

boxplot(Zs_sweep_netting ~ Treatment, data=Cica1,
 main="Z.scutellaris, sweep netting", ylim=ylim)

Catches of Planthoppers and Leafhoppers

Description

Data of a field trial concerning the impact of a genetically modified variety on the abundance of Planthoppers and Leafhoppers. The trial was designed as a randomized complete block design with 8 blocks (Row). In each block, three treatments were randomized: a conventional variety treated with insecticides (Insecticide), a genetically modified variety (GM), and the near-isogenic line (Iso) the to genetically modified line. These data originate from the second year of the trial in Cica1.

Usage

data(Cica2)

Format

A data frame with 24 observations on the following 8 variables.

Field

a factor with 2 levels, 1 2, separating the two major sites of the trial. On field 1, the blocks 1-5 were situated, on field 2, blocks 6-8 were situated.

Row

a factor with 8 levels, specifying the blocks: R1 R2 R3 R4 R5 R6 R7 R8

Year

a numeric vector, with value 2 for the second year

Treatment

a factor with 3 levels, specifying the genetically modified variety GM, the conventional variety treated with insecticides Insecticide, and the variety that was near-isogenic to the GM variety Iso

Au_Bonitur

counts of Auchenorryhncha by visual assessment

Zs_sweep_netting

counts of the major species Zyginidia scutellaris, catched by sweep nets

Zs_yellow_traps

counts of Zyginidia scutellaris, catched by yellow traps

Zs_stick_traps

counts of Zyginidia scutellaris, catched by sticky traps

Source

...

References

See Cica1 for data of the same trial a year earlier

Examples

data(Cica2)


# A comparison of the treatments:

layout(matrix(1:4,ncol=4))

ylim<-range(Cica2[,c("Au_Bonitur","Zs_sweep_netting", "Zs_yellow_traps", "Zs_stick_traps")])

boxplot(Au_Bonitur ~ Treatment, data=Cica2,
 main= "Aucherrhyncha, visual assessment", ylim=ylim, horizontal=TRUE, las=1)

boxplot(Zs_sweep_netting ~ Treatment, data=Cica2,
 main="Z.scutellaris, sweep netting", ylim=ylim, horizontal=TRUE, las=1)

boxplot(Zs_yellow_traps ~ Treatment, data=Cica2,
 main="Z.scutellaris, yellow traps", ylim=ylim, horizontal=TRUE, las=1)

boxplot(Zs_stick_traps ~ Treatment, data=Cica2,
 main="Z.scutellaris, sticky traps", ylim=ylim, horizontal=TRUE, las=1)

# A comparison of sampling methods:

pairs(Cica2[,c("Au_Bonitur","Zs_sweep_netting", "Zs_yellow_traps", "Zs_stick_traps")])

Simulated count data incl. repeated measurements

Description

Simulated data set of repeated count data within subjects.

Usage

data(CountRep)

Format

A data frame with 160 observations on the following 4 variables.

Abundance

a numeric vector with counts simulated from overdispersed and autocorrelated Poisson distributions

ID

a factor with levels N1 N2 ,..., n40, specifying the subject

Time

a factor with levels T1 T2 T3 T4, specifying the time

Treatment

a factor with levels N S1 S2 S3

Examples

data(CountRep)

A simulated data set

Description

A simulated data set, resembling an experiment on the decomposition of plant material from four different varieties in soil.

Usage

data(Decomp)

Format

A data frame with 1152 observations on the following 5 variables.

Block

a numeric vector, giving the names of the Blocks

PID

a numeric vector, giving the the number of the experimental unit

Variety

a factor with levels N S1 S2 S3, specifying the variety assigned to the experimental unit, randomized within each Block

Time

a factor with levels t1 t2 t3 t4 t5 t6, specifying time points at which measurements were taken

Perc

a numeric vector, giving the percentage of material

Details

The experiment is a randomized complete block design, with repeated measurements within each experimental unit and additional subsamplings within each time point. Plant material from four different varieties was deposited in bags in soil of 32 experimental units (coded by the variable PID), where the varieties had been grown in the vegetation period before. A total number of 36 bags was placed in each experimental unit. At six different time points, plant material was excavated and the content of the bags was analysed concerning the percentage of decomposition relative to the content at the begin of the experiment. At each time point, six bags were excavated at each experimental unit. Some bags could not be found anymore (data missing). The objective was to assess whether properties of the varieties obstruct the decompostion of plant material in the soil. The variety N is of special interest, while the varieties S1, S2, S3 are standard varieties.

Note, that this data set merely serves as an example to evaluate clustered data. Neither in the mean effects nor in the actual data points it does resemble a true experiment!

Source

Simulated.

References

Resembling the experimental structure of an experiment performed by Sabine Prescher (JKI Braunschweig).

Examples

data(Decomp)

Soil eklektor data for some families of Diptera

Description

Hatching of some families of Diptera was recorded in summer 2005 using eklektors covering 2 square meters of soil surface each. A total of 32 eklektors were arranged in a randomized field trial. Total counts of individuals over the whole season are reported. Aim was to assess the impact of a novel treatment on the abundance of Diptera with larval development in the soil, compared to three standard treatments.

Usage

data(Diptera)

Format

A data frame with 32 observations on the following 7 variables.

Callip

a numeric vector

Chloro

a numeric vector

Ephyd

a numeric vector

Droso

a numeric vector

Ges

a numeric vector, total number of species

Chiro

a numeric vector

Treatment

a factor, specifying the four different treatments, with levels S1 S2 for two standard treatments, SNovum for the standard treatment most similar to the novel treatment, and Novum, for the novel treatment

Source

personal communications S. Prescher, JKI Braunschweig, Germany

Examples

data(Diptera)


layout(matrix(1:6, nrow=3))

boxplot(Callip~Treatment, data=Diptera, horizontal=TRUE, las=1,
 main="Abundanz Callip", col=c("white","white","blue","red"))

boxplot(Chloro~Treatment, data=Diptera, horizontal=TRUE, las=1,
 main="Abundanz Chloro", col=c("white","white","blue","red"))

boxplot(Ephyd~Treatment, data=Diptera, horizontal=TRUE, las=1,
 main="Abundanz Ephyd", col=c("white","white","blue","red"))

boxplot(Droso~Treatment, data=Diptera, horizontal=TRUE, las=1,
 main="Abundanz Droso", col=c("white","white","blue","red"))

boxplot(Chiro~Treatment, data=Diptera, horizontal=TRUE, las=1,
 main="Abundanz Chiro", col=c("white","white","blue","red"))

boxplot(Ges~Treatment, data=Diptera, horizontal=TRUE, las=1,
 main="Abundanz all Diptera", col=c("white","white","blue","red"))

Simulated example data, drawn from a Negative Binomial Distribution

Description

Simulated example data, response drawn from a Negative Binomial Distribution, Covariables follow a multivariate normal distribution.

Usage

data(ExNBCov)

Format

A data frame with 32 observations on the following 12 variables.

Resp

a numeric vector, a response of counts

Group

a factor with levels A1 A2 A3 A4, e.g. varieties

X1

a numeric covariable

X2

a numeric covariable

X3

a numeric covariable

X4

a numeric covariable

X5

a numeric covariable

X6

a numeric covariable

X7

a numeric covariable

X8

a numeric covariable

X9

a numeric covariable

X10

a numeric covariable

Examples

data(ExNBCov)

boxplot(Resp ~ Group, data=ExNBCov)

pairs(ExNBCov)

Simulated example data following a Poisson distribution

Description

Reponse follows a Poisson distribution, the covariables follow a multivariate normal on the log-scale.

Usage

data(ExPCov)

Format

A data frame with 32 observations on the following 12 variables.

resp

a response of counts

A

a factor with levels A1 A2 A3 A4, e.g. varieties

C1

a numeric covariable

C2

a numeric covariable

C3

a numeric covariable

C4

a numeric covariable

C5

a numeric covariable

C6

a numeric covariable

C7

a numeric covariable

C8

a numeric covariable

C9

a numeric covariable

C10

a numeric covariable

Examples

data(ExPCov)

boxplot(resp ~ A, data=ExPCov)

pairs(ExPCov)

Pupation and Hatching rate in a feeding experiment with four varieties

Description

Larvae of a non-target organism were fed with plant material derived from a novel variety(Novum), material from three standard varieties (NStandard: the standard variety most similar to Novum, and two additional standard varieties S1 and S2). Objective was to assess the impact of Novum on the pupation and hatching rate of an animal that potentially feeds on plant material compared to accepted standard varieties.

Usage

data(Feeding)

Format

A data frame with 32 observations on the following 5 variables.

Rep

a factor with 32 levels indexing the 32 replications

Variety

a factor with 4 levels: S1 and S2 are two standard varieties, Novum is a novel variety, and NStandard is the standard variety most similar to Novum

Total

the total number of animals in each experimental unit

Pupating

number of individuals pupating in each unit, the others died

Hatching

number of individuals hatching from the pupae

Examples

data(Feeding)

# Larval mortality:

Feeding$Lmort <- Feeding$Total - Feeding$Pupating

fit1<-glm(cbind(Pupating,Lmort)~Variety,data=Feeding, family=quasibinomial)
anova(fit1, test="F")

Interaction contrasts for 2-factorial designs

Description

Builds a family of intercation contrasts for complete two-factorial designs.

Usage

IAcontrasts(type, k)

Arguments

Details

Creates contrast matrices using contrMat from package multcomp, creates the kronecker product of both and creates suitable columnnames.

Value

A matrix with k[1]*k[2] columns and a number of rows depending on type.

See Also

for 2-way interaction contrasts directly based on 2 contrasts matrices, see IAcontrastsCMAT; two possibilities to specify appropriate rownames are implemented in function c2compnames

Examples

IAC<-IAcontrasts(type=c("Tukey", "Tukey"), k=c(3,4))

IAC

IAC2<-c2compnames(IAC, ntype="sequ")
IAC2

Interaction contrasts for a two-factorial design

Description

Builds a family of intercation contrasts for complete two-factorial designs.

Usage

IAcontrastsCMAT(CMAT1, CMAT2)

Arguments

Details

Builds the kronecker product of CMAT1 and CMAT2 and creates suitable columnnames. Note that CMAt1 and CMAT2 are not checked, and hence its up to the user to define them suitably.

Value

A matrix with k[1]*k[2] columns.

See Also

for interaction contrasts based on contrast definition and the number of levels of the factors in atwo-way layout, see IAcontrasts; two possibilities to specify appropriate rownames are implemented in function c2compnames

Examples

library(multcomp)


n1<-c(10,10,10,10)
names(n1)<-c("A","B","C","D")

n2<-c(3,3,3)
names(n2)<-c(1,2,3)

CMT1<-contrMat(n1, type="Tukey")

CMT2<-contrMat(n2, type="Tukey")

IAC<-IAcontrastsCMAT(CMAT1=CMT1, CMAT2=CMT2)

c2compnames(IAC, ntype="sequ")

###


n1<-c(10,10,10,10)
names(n1)<-c("A","B","C","D")

n2<-c(3,3,3)
names(n2)<-c(1,2,3)

CMD1<-contrMat(n1, type="Dunnett")

CMD2<-contrMat(n2, type="Dunnett")

IAC<-IAcontrastsCMAT(CMAT1=CMD1, CMAT2=CMD2)

c2compnames(IAC, ntype="sequ")

Insect counts of 12 Species

Description

Simulated data, inpired by a real field investigating the potential impact of genetically modified crop on several insect species belonging to the same order. The trial was designed as a randomized complete block design with 8 blocks (Block), and a total of 24 plots. In each block, three treatments (Treatment) were randomized: a conventional variety treated with insecticides (Ins), a genetically modified variety (GM) without insecticide treatment, and the near-isogenic variety (Iso) the to genetically modified variety, without insecticide treatment. Individuals were counted (after classification to the species level) in two different dates in each year of the trial, where the the second date was of higher importance for assessment of impacts of GM variety on non-target species. In total 12 Species were observed during the trial.

Usage

data(Lepi)

Format

A data frame with 144 observations on the following 17 variables.

Year

a numeric vector, the year 1, 2, 3

Date

a numeric vector, 1 and 2 separating the 2 sampling date in each year

Block

a numeric vector, with values 1-8, indicator variable for the 8 blocks

Treatment

a factor with three levels identifying the varieties: GM is the genetically modified variety, Ins the conventional variety with insecticide treatment and Iso the near isognic line without insecticide treatment

Plot

a factor with 24 levels, identifying the individual plots

Sp1

counts of taxon 1

Sp2

counts of taxon 2

Sp3

counts of taxon 3

Sp4

counts of taxon 4

Sp5

counts of taxon 5

Sp6

counts of taxon 6

Sp7

counts of taxon 7

Sp8

counts of taxon 8

Sp9

counts of taxon 9

Sp10

counts of taxon 10

Sp11

counts of taxon 11

Sp12

counts of taxon 12

Source

Simulated data.

Examples

data(Lepi)

str(Lepi)

summary(Lepi)

SPEC<-names(Lepi)[-(1:5)]

# Occurrence

occur<-lapply(X=Lepi[,SPEC], FUN=function(x){length(which(x>0))})

unlist(occur)

# Species with reasonable occurence in the whole data:

SPEC2<-SPEC[c(1,2,3,6,8,9,11)]

pairs(Lepi[,SPEC2])

# 


layout(matrix(1:2, ncol=1 ))
par(mar=c(2,8,2,1))

boxplot(Sp2 ~ Treatment*Year, data=Lepi, main="Species 2",
 las=1, horizontal=TRUE, col=c("red","white","white"))

boxplot(Sp3 ~ Treatment*Year, data=Lepi, main="Species 3",
 las=1, horizontal=TRUE, col=c("red","white","white"))

Simulated data set for a simple mixed model

Description

Simulated data set for a simple mixed model

Usage

data(MM1)

Format

A data frame with 160 observations on the following 3 variables.

Y

a numeric vector, the response, sampled from a normal distribution

F

a factor with levels F1 F2 F3 F4, representing fixed effects

R

a factor with levels R1 R2 R3 R4 R5, representing random effects, sampled from a normal distribution

Examples

data(MM1)
boxplot(Y~F*R, data=MM1)

Simulated data for a simple mixed model with Poisson response

Description

A fixed factor F (four levels) and a random factor (five levels), modifying the mean response (random Intercept) Y is a variable following a Poisson distribution.

Usage

data(MMPois)

Format

A data frame with 160 observations on the following 3 variables.

Y

a numeric vector, the Poisson distributed response.

F

a factor with levels F1 F2 F3 F4

R

a factor with levels R1 R2 R3 R4 R5

Source

simulation

Examples

data(MMPois)

boxplot(Y~R*F, data=MMPois, las=2)

Simulated data for a simple mixed model with Poisson response

Description

Simulated data with a fixed factor cult (4 levels), with 8 randomized replications each, a (fixed) factor time (6 levels), which are repeated measurements taken from the same experimental units. The 32 experimental (plotid) units differ in their mean response follwing a gaussian distribution. The response Y follows a Poisson distribution.

Usage

data(MMPoisRep)

Format

A data frame with 192 observations on the following 4 variables.

plotid

a factor with 32 levels, representing the 32 experimental units (plots)

cult

a factor with 4 levels (C1 C2 C3 C4), representing a fixed factor (e.g. the cultivar)

time

a factor with 6 levels (T1 T2 T3 T4 T5 T6) specifying repeated measurements on the same experimental units (plotid) over time

Y

a numeric vector, following a Poisson distribution

Source

simulation

Examples

data(MMPoisRep)

boxplot(Y ~ cult*time, data=MMPoisRep, las=TRUE)

Trap counts of Nematocera

Description

In a field trial, 4 treatments (A,B, Standard, Novum) were arranged in a randomized complete block design with 8 blocks and 32 plots. In summer 2005 soil eklektor traps were placed in each plot, on six dates from 2005-07-12 to 2005-09-25, the number of individuals of Nematocera (gnats, midges and others) hatching from soil were counted. The individuals were classified to the family level. Interest was in assessing potential effects of a novel agricultural practice (Novum) on the abundance of Nematocera.

Usage

data(Nematocera)

Format

A data frame with 192 observations on the following 14 variables.

Date

a POSIXt, the time of counting the individuals in the eklektor trap

Treatment

a factor with 4 levels, A, B, Standard and Novum, where Novum is the novel treatment, Standard is the standard treatment most similar to Novum, and A and B are additional standard treatments.

Block

a numeric vector, specifying the blocks 1-8

Plot

a factor with 32 levels A1 to Standard8, indicator variables for the individual eklektors

Bibio

a numeric vector, counts of individuals, belonging to the family

Cecido

a numeric vector, counts of individuals, belonging to the family

Cerato

a numeric vector, counts of individuals, belonging to the family

Chiro

a numeric vector, counts of individuals, belonging to the family

Myceto

a numeric vector, counts of individuals, belonging to the family

Psycho

a numeric vector, counts of individuals, belonging to the family

Scato

a numeric vector, counts of individuals, belonging to the family

Sciari

a numeric vector, counts of individuals, belonging to the family

Tipuli

a numeric vector, counts of individuals, belonging to the family

Total

a numeric vector, total count of individuals belonging to the suborder Nematocera

Source

personal communications, S.Prescher, JKI Braunschweig, Germany

Examples

data(Nematocera)

par(mar=c(11,5,3,1))

boxplot(Total ~ Treatment*Date, data=Nematocera, las=2, col=c("white","white","blue","green"))
legend(x=15, y=100, legend=levels(Nematocera$Treatment), fill=c("white","white","blue","green"))

pairs(Nematocera[,c("Cecido","Cerato","Chiro","Myceto","Psycho","Sciari")])

For internal use

Description

Transform a data set to a dataset appropriate for certain OpenBUGS models.

Usage

R2Bugsdat1w(formula, data)

Arguments

Details

For internal use.

Value

a list, containing the elements

For internal use.

Description

Transform a data set to a dataset appropriate for certain OpenBUGS models.

Usage

R2Bugsdat1w.data.frame(data, response, treatment, Intercept = FALSE)

Arguments

Details

For internal use.

Value

a list, containing the elements

Simultaneous confidence sets from empirical joint distribution.

Description

Calcualte simultaneous confidence sets according to Besag et al. (1995) from a empirical joint distribution of a parameter vector. Joint empirical distributions might be obtained from WinBUGS or OpenBUGS calls.

Usage

## Default S3 method:
SCSnp(x, conf.level = 0.95,
 alternative = "two.sided", ...)

## S3 method for class 'bugs'
SCSnp(x, conf.level = 0.95,
 alternative = "two.sided", whichp = NULL, ...)

## S3 method for class 'CCRatio'
SCSnp(x, ...)

## S3 method for class 'CCDiff'
SCSnp(x, ...)

Arguments

Details

Let P be the number of parameters in the parameter vector and N be the total number of values obtained for the empirical joint distribution of the parameter vector, e.g. as can be obtaine e.g., from Gibbs sampling.

Value

An object of class "SCSnp", a list with elements

Author(s)

Frank Schaarschmidt, adapting an earliere version of Gemechis D. Djira

References

Besag J, Green P, Higdon D, Mengersen K (1995): Bayesian Computation and Stochastic Systems. Statistical Science 10 (1), 3-66.

See Also

CInp for a wrapper to quantile to compute simple percentile intervals on each of P marginal distributions

Examples

# Assume a 1000 times 4 matrix of 4 mutually independent
# normal variables:

X<-cbind(rnorm(1000), rnorm(1000), rnorm(1000), rnorm(1000))

SCSts<-SCSnp(x=X, conf.level=0.9, alternative="two.sided")
SCSts

SCS<-SCSts$conf.int

in1<-X[,1]>=SCS[1,1] & X[,1]<=SCS[1,2] 

in2<-X[,2]>=SCS[2,1] & X[,2]<=SCS[2,2] 

in3<-X[,3]>=SCS[3,1] & X[,3]<=SCS[3,2] 

in4<-X[,4]>=SCS[4,1] & X[,4]<=SCS[4,2] 

sum(in1*in2*in3*in4)

Transform confidence intervals from glm fits.

Description

Transform confidence intervals derived from glm fits back to original scale and give appropriate names.

Usage

## S3 method for class 'glht'
UnlogCI(x)

Arguments

Details

Applies exponential function on the estimates and confidence limits and creates useful names for the comparisons and parameters.

Value

An object of class "UnlogCI".

See Also

plotCI.UnlogCI for plotting the result

Examples

# # # CI for odds ratios
# # # for models on the logit-link

data(Feeding)

# Larval mortality:

Feeding$Lmort <- Feeding$Total - Feeding$Pupating

fit1<-glm(cbind(Pupating,Lmort)~Variety,data=Feeding, family=quasibinomial)
anova(fit1, test="F")

library(multcomp)

comp<-glht(fit1, mcp(Variety="Tukey"))

CIraw<-CIGLM(comp,method="Raw")

CIraw

UnlogCI(CIraw)

plotCI(UnlogCI(CIraw), lines=c(0.25,0.5,2,4),
 lineslwd=c(1,2,2,1), linescol=c("red","black","black","red"))


# # # # # # #

# # #  CI for ratios of means
# # # for models on the log-link

data(Diptera)

# Larval mortality:

fit2<-glm(Ges~Treatment, data=Diptera, family=quasipoisson)
anova(fit2, test="F")

library(multcomp)

comp<-glht(fit2, mcp(Treatment="Tukey"))

CIadj<-CIGLM(comp,method="Adj")

CIadj

UnlogCI(CIadj)

plotCI(UnlogCI(CIadj), lines=c(0.5,1,2), lineslwd=c(2,1,1))

Allignment according to one factor

Description

Substracts the mean or median from observations belonging to the same level of a factor.

Usage

allignment(response, block, type = c("mean", "median"), ...)

Arguments

Details

Splits response according to the levels of block, calculates and substracts the mean or median and returns the resulting vector in appropriate order.

Value

A numeric vector.

Define row names of a contrast matrix, depending on its column names

Description

Define row names of a contrast matrix, depending on its column names, as can be necessary for contrasts matrices. Currently, two options to do that are available.

Usage

c2compnames(cmat, ntype = "aggr")

Arguments

Value

The input matrix cmat, with its row names replaced.

See Also

contrMat in multcomp to define contrast matrices of different types

Examples

# names for interaction contrasts:

n1<-c(10,10,10,10)
names(n1)<-c("A","B","C","D")

n2<-c(3,3,3)
names(n2)<-c(1,2,3)

library(multcomp)

CMT1<-contrMat(n1, type="Tukey")

CMT2<-contrMat(n2, type="Tukey")

IAC<-IAcontrastsCMAT(CMAT1=CMT1, CMAT2=CMT2)

c2compnames(IAC, ntype="aggr")

c2compnames(IAC, ntype="sequ")


###############################

# names for Williams-type contrasts:

n1<-c(10,10,10,10)
names(n1)<-c("C0","D1","D5","D10")


CMW<-contrMat(n1, type="Williams")

CMW

c2compnames(CMW, ntype="aggr")

c2compnames(CMW, ntype="sequ")

For internal use.

Description

For internal use. Check the appropriateness of input arguments of function simplesimint.

Usage

checkargssim(coef, vcov, cmat)

Arguments

A simulated data set of lognormal data

Description

A simulated data set of lognormal data, could be concentrations

Usage

data(fakeln)

Format

A data frame with 32 observations on the following 2 variables.

Concmug

a numeric vector, serving as response variable

Treat

a factor with levels N S Sa Sb

Examples

data(fakeln)
boxplot(Concmug ~ Treat, data=fakeln)

Replace - by / in character strings

Description

Replace - by / in character strings

Usage

minus2slash(x)

Arguments

Value

a vector of character strings

For internal use. Extract model parameters needed for multcomp from objects of class gamlss or geeglm.

Description

Only for internal use with glht in package multcomp. Extracts model parameters needed for glht(multcomp) from objects of class gamlss or geeglm.

Usage

## S3 method for class 'gamlss'
modelparm(model, coef.=coef, vcov.=vcov, df = NULL, ...)
## S3 method for class 'geeglm'
modelparm(model, coef.=coef, vcov.=vcov, df = NULL, ...)

Arguments

Details

Test versions. Only for internal use. Further checks still need to be implemented.

Value

As modelparm.default in package multcomp, a list with elements

Author(s)

Daniel Gerhard, Frank Schaarschmidt adapting code by T.Hothorn available in modelparm.default, package multcomp

See Also

function glht in package multcomp

Plot confidence intervals calculated by pairwiseCI

Description

Plot confidence intervals calculated by calling pairwiseCI or UnlogCI.

Usage

## S3 method for class 'UnlogCI'
plotCI(x, ...)
## S3 method for class 'simplesimint'
plotCI(x, ...)

Arguments

Value

A plot.

Print function for SCSnp

Description

A print functions.

Usage

## S3 method for class 'SCSnp'
print(x, ...)
## S3 method for class 'CInp'
print(x, ...)
## S3 method for class 'simplesimint'
print(x, ...)
## S3 method for class 'UnlogCI'
print(x, ...)

Arguments

Simultaneous confidence intervals from raw estimates

Description

Calculates simultaneous confidence intervals for multiple contrasts based on a parameter vector, its variance-covariance matrix and (optionally) the degrees of freedom, using quantiles of the multivar

Usage

simplesimint(coef, vcov, cmat, df = NULL, conf.level = 0.95,
 alternative = c("two.sided", "less", "greater"))

Arguments

Details

Implements the methods formerly available in package multcomp, function csimint. Input values are a vector of parameter estimates \mu of length P, a corresponding estimate for its variance-covariance matrix \Sigma (P times P), and a contrast matrix C of dimension M \times P. The contrasts L = C \mu are computed, the variance-covariance matrix (being a function of C and \Sigma) and the corresponding correlation matrix R are computed. Finally, confidence intervals for L are computed: if df is given, quantiles of an M-dimensional t distribution with correlation matrix R are used, otherwise quantiles of an M-dimensional standard normal distribution with correlation matrix R are used.

Value

An object of class "simplesimint"

Note

This is a testversion and has not been checked extensively.

Author(s)

Frank Schaarschmidt

See Also

See ?coef and ?vcov for extracting of parameter vectors and corresponding variance covariance matrices from various model fits.

Examples

# For the simple case of Gaussian response
# variables with homoscedastic variance,
# see the following example


library(mratios)
data(angina)

boxplot(response ~ dose, data=angina)

# Fit a cell means model,

fit<-lm(response ~ 0+dose, data=angina)

# extract cell means, the corresponding
# variance-covariance matrix and the
# residual degree of freedom,

cofi<-coef(fit)
vcofi<-vcov(fit)
dofi<-fit$df.residual

# define an appropriate contrast matrix,
# here, comparisons to control

n<-unlist(lapply(split(angina$response, f=angina$dose), length))
names(n)<-names(cofi)

cmat<-contrMat(n=n, type="Dunnett")
cmat

#

test<-simplesimint(coef=cofi, vcov=vcofi, df=dofi, cmat=cmat, alternative="greater" )

test

summary(test)

plotCI(test)

### Note, that the same result can be achieved much more conveniently
### using confint.glht in package multcomp

Detailed print out for simplesimint objects

Description

Produce a detailed print out for objects of class "simplesimint".

Usage

## S3 method for class 'simplesimint'
summary(object, ...)

Arguments

Value

Print

Extract variance covariance matrix from objects of class gamlss

Description

Only for internal use. Extract the covariance matrix corresponding to the mu parameters of a gamlss-fit.

Usage

## S3 method for class 'gamlss'
vcov(object, ...)

Arguments

Details

Only for internal use. Needs implementation of warnings.

Value

A matrix of dimension mxm, if m is the length of mu-parameters from a gamlss fit.

Author(s)

Daniel Gerhard

See Also

packages gamlss and multcomp

Extract variance covariance matrix from objects of class geeglm

Description

Only for internal use with function glht: Extract the variance covariance matrix corresponding to the mu parameters of a gamlss-fit. Merely a wrapper of the method internally used in function summary.geeglm, package geepack.

Usage

## S3 method for class 'geeglm'
vcov(object, ...)

Arguments

Details

Test version. Only for internal use. Needs implementation of warnings.

Value

A matrix of dimension m times m, if m is the length of coefficients from a geeglm fit.

See Also

packages gamlss and multcomp