Title: | Safety Assessment in Agricultural Field Trials |
---|---|
Description: | Collection of functions, data sets and code examples for evaluations of field trials with the objective of equivalence assessment. |
Authors: | Frank Schaarschmidt |
Maintainer: | Frank Schaarschmidt <[email protected]> |
License: | GPL-2 |
Version: | 0.1-10 |
Built: | 2024-11-19 06:28:42 UTC |
Source: | CRAN |
Substracts the mean or median from observations belonging to the same level of a factor.
allignment(response, block, type = c("mean", "median"), ...)
allignment(response, block, type = c("mean", "median"), ...)
response |
a numeric vector |
block |
a factor of the same length as |
type |
type of location measure to calculate and substract; only the choices "mean" and "median" are supported |
... |
further arguments to be passed to |
Splits response
according to the levels of block
, calculates and substracts the mean or median and returns the resulting vector in appropriate order.
A numeric vector.
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.
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, ...)
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, ...)
X |
a data.frame of dimension n times p, containing integer entries as species counts of p species from n independent samplings |
f |
a factor, usually with more than two levels. Must be of length n, when X is an n times p matrix |
type |
a single character string, defining a contrast type. Supported options are 'Dunnett','Tukey','Sequen'; for more options, see |
cmat |
user defined contrasts:
when using |
conf.level |
a single numeric value between 0.5 and 1 |
alternative |
a single character string, one of 'two.sided','less' and 'greater' |
madj |
a single logical value, indicating whether simultaneous (if |
... |
Further arguments to be passed to the function |
NOTE: This is a test version!
If madj=TRUE
, an object of class "SCSnp", see SCSnp
for details.
If madj=FALSE
, an object of class "CInp", see CInp
for details.
Frank Schaarschmidt
SCSnp
,
these function internally make use of CCDiff.boot
, CCDiff.default
,
CCRatio.boot
, CCRatio.default
,
boot
and estSimpsonf
.
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
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
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.
data(Brachycera)
data(Brachycera)
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
...
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"))
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"))
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.
c2compnames(cmat, ntype = "aggr")
c2compnames(cmat, ntype = "aggr")
cmat |
a contrast matrix |
ntype |
a single character string,
defining how to build names from the column names of |
The input matrix cmat, with its row names replaced.
contrMat
in multcomp to define contrast matrices of different types
# 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")
# 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")
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.
data(Cica1)
data(Cica1)
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
...
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)
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)
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.
data(Cica2)
data(Cica2)
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
...
See Cica1
for data of the same trial a year earlier
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")])
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")])
Computes confidence intervals from the output of a glm, by calling to glht(multcomp).
CIGLM(x, conf.level = 0.95, method = c("Raw", "Adj", "Bonf"))
CIGLM(x, conf.level = 0.95, method = c("Raw", "Adj", "Bonf"))
x |
a object of class |
conf.level |
confidence level, a single numeric value between 0.5 and 1 |
method |
a single character string, with |
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!
An object of class "confint.glht"
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.
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
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 for each parameter from the empirical joint distribution of a parameter vector of length P.
## 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, ...)
## 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, ...)
x |
an N-times-P matrix, or an object of class |
conf.level |
a single numeric value between 0.5 and 1, specifying the local confidence level for each of the P parameters |
alternative |
a single character string, one of |
whichp |
a single character string, naming an element of the |
... |
currently not used |
Construct simple confidence intervals based on order statistics applied to the marginal empirical distributions in x
.
An object of class "CInp", a list with elements
conf.int |
a P-times-2 matrix containing the lower and upper confidence limits |
estimate |
a numeric vector of length P, containing the medians of the P marginal empirical distributions |
x |
the input object |
k |
the number of values outside each confidence interval, i.e. conf.level*N |
N |
the number of values used to construct each confidence interval |
conf.level |
a single numeric value, the nominal confidence level, as input |
alternative |
a single character string, as input |
The function internally used is quantile
with its default settings.
See SCSnp
for simultaneous sets.
# 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] ) )
# 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] ) )
Simulated data set of repeated count data within subjects.
data(CountRep)
data(CountRep)
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
data(CountRep)
data(CountRep)
A simulated data set, resembling an experiment on the decomposition of plant material from four different varieties in soil.
data(Decomp)
data(Decomp)
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
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!
Simulated.
Resembling the experimental structure of an experiment performed by Sabine Prescher (JKI Braunschweig).
data(Decomp)
data(Decomp)
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.
data(Diptera)
data(Diptera)
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
personal communications S. Prescher, JKI Braunschweig, Germany
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"))
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, response drawn from a Negative Binomial Distribution, Covariables follow a multivariate normal distribution.
data(ExNBCov)
data(ExNBCov)
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
data(ExNBCov) boxplot(Resp ~ Group, data=ExNBCov) pairs(ExNBCov)
data(ExNBCov) boxplot(Resp ~ Group, data=ExNBCov) pairs(ExNBCov)
Reponse follows a Poisson distribution, the covariables follow a multivariate normal on the log-scale.
data(ExPCov)
data(ExPCov)
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
data(ExPCov) boxplot(resp ~ A, data=ExPCov) pairs(ExPCov)
data(ExPCov) boxplot(resp ~ A, data=ExPCov) pairs(ExPCov)
A simulated data set of lognormal data, could be concentrations
data(fakeln)
data(fakeln)
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
data(fakeln) boxplot(Concmug ~ Treat, data=fakeln)
data(fakeln) boxplot(Concmug ~ Treat, data=fakeln)
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.
data(Feeding)
data(Feeding)
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
data(Feeding) # Larval mortality: Feeding$Lmort <- Feeding$Total - Feeding$Pupating fit1<-glm(cbind(Pupating,Lmort)~Variety,data=Feeding, family=quasibinomial) anova(fit1, test="F")
data(Feeding) # Larval mortality: Feeding$Lmort <- Feeding$Total - Feeding$Pupating fit1<-glm(cbind(Pupating,Lmort)~Variety,data=Feeding, family=quasibinomial) anova(fit1, test="F")
Builds a family of intercation contrasts for complete two-factorial designs.
IAcontrasts(type, k)
IAcontrasts(type, k)
type |
a vector of two character strings, specifying the contrast |
k |
a vector of two integers, specifying the number of groups in each factor of the two-factorial design |
Creates contrast matrices using contrMat
from package multcomp, creates the kronecker product of both and creates suitable columnnames.
A matrix with k[1]*k[2] columns and a number of rows depending on type.
for 2-way interaction contrasts directly based on 2 contrasts matrices, see IAcontrastsCMAT
;
two possibilities to specify appropriate rownames are implemented in function c2compnames
IAC<-IAcontrasts(type=c("Tukey", "Tukey"), k=c(3,4)) IAC IAC2<-c2compnames(IAC, ntype="sequ") IAC2
IAC<-IAcontrasts(type=c("Tukey", "Tukey"), k=c(3,4)) IAC IAC2<-c2compnames(IAC, ntype="sequ") IAC2
Builds a family of intercation contrasts for complete two-factorial designs.
IAcontrastsCMAT(CMAT1, CMAT2)
IAcontrastsCMAT(CMAT1, CMAT2)
CMAT1 |
a (named) contrast matrix |
CMAT2 |
a (named) contrast matrix |
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.
A matrix with k[1]*k[2] columns.
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
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")
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")
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.
data(Lepi)
data(Lepi)
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
Simulated data.
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"))
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
data(MM1)
data(MM1)
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
data(MM1) boxplot(Y~F*R, data=MM1)
data(MM1) boxplot(Y~F*R, data=MM1)
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.
data(MMPois)
data(MMPois)
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
simulation
data(MMPois) boxplot(Y~R*F, data=MMPois, las=2)
data(MMPois) boxplot(Y~R*F, data=MMPois, las=2)
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.
data(MMPoisRep)
data(MMPoisRep)
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
simulation
data(MMPoisRep) boxplot(Y ~ cult*time, data=MMPoisRep, las=TRUE)
data(MMPoisRep) boxplot(Y ~ cult*time, data=MMPoisRep, las=TRUE)
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.
data(Nematocera)
data(Nematocera)
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
personal communications, S.Prescher, JKI Braunschweig, Germany
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")])
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")])
Plot confidence intervals calculated by calling pairwiseCI or UnlogCI.
## S3 method for class 'UnlogCI' plotCI(x, ...) ## S3 method for class 'simplesimint' plotCI(x, ...)
## S3 method for class 'UnlogCI' plotCI(x, ...) ## S3 method for class 'simplesimint' plotCI(x, ...)
x |
an object of class |
... |
further arguments to be passed to |
A plot.
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.
## 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, ...)
## 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, ...)
x |
a matrix N-times-P matrix or an object of class |
conf.level |
a single numeric value between 0.5 and 1, the simultaneous confidence level |
alternative |
a single character string, one of |
whichp |
a single character string, naming an element of the |
... |
further arguments, currently not used |
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.
An object of class "SCSnp", a list with elements
conf.int |
a P-times-2 matrix containing the lower and upper confidence limits |
estimate |
a numeric vector of length P, containing the medians of the P marginal empirical distributions |
x |
the input object |
k |
the number of values outside the SCS, i.e. conf.level*N |
N |
the number of values used to construct the confidence set |
conf.level |
a single numeric value, the nominal confidence level, as input |
alternative |
a single character string, as input |
Frank Schaarschmidt, adapting an earliere version of Gemechis D. Djira
Besag J, Green P, Higdon D, Mengersen K (1995): Bayesian Computation and Stochastic Systems. Statistical Science 10 (1), 3-66.
CInp
for a wrapper to quantile
to compute simple percentile intervals on each of P marginal distributions
# 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)
# 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)
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
simplesimint(coef, vcov, cmat, df = NULL, conf.level = 0.95, alternative = c("two.sided", "less", "greater"))
simplesimint(coef, vcov, cmat, df = NULL, conf.level = 0.95, alternative = c("two.sided", "less", "greater"))
coef |
a single numeric vector, specifying the point estimates of the parameters of interest |
vcov |
the variance-covariance matrix corresponding to |
cmat |
the contrasts matrix specifying the comparisons of interest with respect to |
df |
optional, the degree of freedom for the multivariate t-distribution; if specified, quantiles from the multivariate t-distribution are used for confidence interval estimation, if not specified (default), quantiles of the multivariate normal distribution are used |
conf.level |
a single numeric value between 0.5 and 1.0; the simultaneous confidence level |
alternative |
a single character string, |
Implements the methods formerly available in package multcomp, function csimint
.
Input values are a vector of parameter estimates of length
,
a corresponding estimate for its variance-covariance matrix
(P times P), and a
contrast matrix
of dimension
. The contrasts
are computed,
the variance-covariance matrix (being a function of
and
) and the corresponding correlation matrix
are computed.
Finally, confidence intervals for
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.
An object of class "simplesimint"
estimate |
the estimates of the contrasts |
lower |
the lower confidence limits |
upper |
the upper confidence limits |
cmat |
the contrast matrix, as input |
alternative |
a character string, as input |
conf.level |
a numeric value, as input |
quantile |
a numeric value, the quantile used for confidence interval estimation |
df |
a numeric value or NULL, as input |
stderr |
the standard error of the contrasts |
vcovC |
the variance covariance matrix of the contrasts |
This is a testversion and has not been checked extensively.
Frank Schaarschmidt
See ?coef
and ?vcov
for extracting of parameter vectors and corresponding variance covariance matrices from various model fits.
# 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
# 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
Produce a detailed print out for objects of class "simplesimint".
## S3 method for class 'simplesimint' summary(object, ...)
## S3 method for class 'simplesimint' summary(object, ...)
object |
an object of class |
... |
further arguments to be passed to |
Transform confidence intervals derived from glm fits back to original scale and give appropriate names.
## S3 method for class 'glht' UnlogCI(x)
## S3 method for class 'glht' UnlogCI(x)
x |
an object of class |
Applies exponential function on the estimates and confidence limits and creates useful names for the comparisons and parameters.
An object of class "UnlogCI"
.
plotCI.UnlogCI
for plotting the result
# # # 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))
# # # 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))
Only for internal use. Extract the covariance matrix corresponding to the mu parameters of a gamlss-fit.
## S3 method for class 'gamlss' vcov(object, ...)
## S3 method for class 'gamlss' vcov(object, ...)
object |
An object of class "gamlss" as can be created by calling |
... |
Currently not used. |
Only for internal use. Needs implementation of warnings.
A matrix of dimension mxm, if m is the length of mu-parameters from a gamlss fit.
Daniel Gerhard
packages gamlss
and multcomp
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.
## S3 method for class 'geeglm' vcov(object, ...)
## S3 method for class 'geeglm' vcov(object, ...)
object |
An object of class "geeglm" as can be created by calling |
... |
Currently not used. |
Test version. Only for internal use. Needs implementation of warnings.
A matrix of dimension m times m, if m is the length of coefficients from a geeglm fit.
packages gamlss
and multcomp