Package 'glmtoolbox'

Description

Computes the adjusted R-squared

Usage

adjR2(..., digits, verbose)

Arguments

Value

A matrix with the values of the adjusted R-squared for all model fit objects.

Description

Computes the adjusted deviance-based R-squared in generalized linear models.

Usage

## S3 method for class 'glm'
adjR2(..., digits = max(3, getOption("digits") - 2), verbose = TRUE)

Arguments

Details

The deviance-based R-squared is computed as $R^2=1 - Deviance/Null.Deviance$. Then, the adjusted deviance-based R-squared is computed as $1 - \frac{n-1}{n-p}(1-R^2)$, where $p$ is the number of parameters in the linear predictor and $n$ is the sample size.

Value

a matrix with the following columns

Examples

###### Example 1: Fuel efficiency of cars
Auto <- ISLR::Auto
fit1 <- glm(mpg ~ horsepower*weight, family=Gamma(inverse), data=Auto)
fit2 <- update(fit1, formula=mpg ~ horsepower*weight*cylinders)
fit3 <- update(fit1, family=Gamma(log))
fit4 <- update(fit2, family=Gamma(log))
fit5 <- update(fit1, family=inverse.gaussian(log))
fit6 <- update(fit2, family=inverse.gaussian(log))

AIC(fit1,fit2,fit3,fit4,fit5,fit6)
BIC(fit1,fit2,fit3,fit4,fit5,fit6)
adjR2(fit1,fit2,fit3,fit4,fit5,fit6)

Description

Computes the adjusted deviance-based R-squared in generalized nonlinear models.

Usage

## S3 method for class 'gnm'
adjR2(..., digits = max(3, getOption("digits") - 2), verbose = TRUE)

Arguments

Details

The deviance-based R-squared is computed as $R^2=1 - Deviance/Null.Deviance$. Then, the adjusted deviance-based R-squared is computed as $1 - \frac{n-1}{n-p}(1-R^2)$, where $p$ is the number of parameters in the "linear" predictor and $n$ is the sample size.

Value

a matrix with the following columns

Examples

###### Example 1: The effects of fertilizers on coastal Bermuda grass
data(Grass)
fit1 <- gnm(Yield ~ b0 + b1/(Nitrogen + a1) + b2/(Phosphorus + a2) + b3/(Potassium + a3),
            family=gaussian(inverse), start=c(b0=0.1,b1=13,b2=1,b3=1,a1=45,a2=15,a3=30), data=Grass)
fit2 <- update(fit1, family=Gamma(inverse))
adjR2(fit1,fit2)

###### Example 2: Developmental rate of Drosophila melanogaster
data(Drosophila)
fit1 <- gnm(Duration ~ b0 + b1*Temp + b2/(Temp-a), family=Gamma(log),
            start=c(b0=3,b1=-0.25,b2=-210,a=55), weights=Size, data=Drosophila)
fit2 <- update(fit1, family=inverse.gaussian(log))
adjR2(fit1,fit2)

###### Example 3: Radioimmunological Assay of Cortisol
data(Cortisol)
fit1 <- gnm(Y ~ b0 + (b1-b0)/(1 + exp(b2+ b3*lDose))^b4, family=Gamma(identity),
            start=c(b0=130,b1=2800,b2=3,b3=3,b4=0.5), data=Cortisol)
fit2 <- update(fit1, family=gaussian(identity))
adjR2(fit1,fit2)

###### Example 4: Age and Eye Lens Weight of Rabbits in Australia
data(rabbits)
fit1 <- gnm(wlens ~ b1 - b2/(age + b3), family=Gamma(log),
            start=c(b1=5.5,b2=130,b3=35), data=rabbits)
fit2 <- update(fit1, family=gaussian(log))
adjR2(fit1,fit2)

Description

Extracts the adjusted R-squared in normal linear models.

Usage

## S3 method for class 'lm'
adjR2(..., digits = max(3, getOption("digits") - 2), verbose = TRUE)

Arguments

Details

The R-squared is computed as $R^2=1 - RSS/Null.RSS$. Then, the adjusted R-squared is computed as $1 - \frac{n-1}{n-p}(1-R^2)$, where $p$ is the number of parameters in the linear predictor and $n$ is the sample size.

Value

a matrix with the following columns

Examples

###### Example 1: Fuel efficiency of cars
fit1 <- lm(mpg ~ log(hp) + log(wt) + qsec, data=mtcars)
fit2 <- lm(mpg ~ log(hp) + log(wt) + qsec + log(hp)*log(wt), data=mtcars)
fit3 <- lm(mpg ~ log(hp)*log(wt)*qsec, data=mtcars)

AIC(fit1,fit2,fit3)
BIC(fit1,fit2,fit3)
adjR2(fit1,fit2,fit3)

Description

The advertising data set consists of sales of that product in 200 different markets. It also includes advertising budgets for the product in each of those markets for three different media: TV, radio, and newspapers.

Usage

data(advertising)

Format

A data frame with 200 rows and 4 variables:

TV

a numeric vector indicating the advertising budget on TV.

radio

a numeric vector indicating the advertising budget on radio.

newspaper

a numeric vector indicating the advertising budget on newspaper.

sales

a numeric vector indicating the sales of the interest product.

Source

https://www.statlearning.com/s/Advertising.csv

References

James G., Witten D., Hastie T., Tibshirani R. (2013, page 15) An Introduction to Statistical Learning with Applications in R, Springer, New York.

Examples

data(advertising)
breaks <- with(advertising,quantile(radio,probs=c(0:3)/3))
labels <- c("low","mid","high")
advertising2 <- within(advertising,radioC <- cut(radio,breaks,labels,include.lowest=TRUE))
dev.new()
with(advertising2,plot(TV,sales,pch=16,col=as.numeric(radioC)))
legend("topleft",legend=c("low","mid","high"),fill=c(1:3),title="Radio",bty="n")

Description

Computes the Akaike-type penalized Gaussian pseudo-likelihood criterion (AGPC) for one or more objects of the class glmgee.

Usage

AGPC(..., k = 2, verbose = TRUE, digits = max(3, getOption("digits") - 2))

Arguments

Details

If k is set to 0 then the AGPC reduces to the Gaussian pseudo-likelihood criterion (GPC), proposed by Carey and Wang (2011), which corresponds to the logarithm of the multivariate normal density function.

Value

A data.frame with the values of the gaussian pseudo-likelihood, the number of parameters in the linear predictor plus the number of parameters in the correlation matrix, and the value of AGPC for each glmgee object in the input.

References

Carey V.J., Wang Y.-G. (2011) Working covariance model selection for generalized estimating equations. Statistics in Medicine 30:3117-3124.

Zhu X., Zhu Z. (2013) Comparison of Criteria to Select Working Correlation Matrix in Generalized Estimating Equations. Chinese Journal of Applied Probability and Statistics 29:515-530.

Fu L., Hao Y., Wang Y.-G. (2018) Working correlation structure selection in generalized estimating equations. Computational Statistics 33:983-996.

Vanegas L.H., Rondon L.M., Paula G.A. (2023) Generalized Estimating Equations using the new R package glmtoolbox. The R Journal 15:105-133.

See Also

QIC, CIC, RJC, GHYC, SGPC

Examples

###### Example 1: Effect of ozone-enriched atmosphere on growth of sitka spruces
data(spruces)
mod1 <- size ~ poly(days,4) + treat
fit1 <- glmgee(mod1, id=tree, family=Gamma(log), data=spruces)
fit2 <- update(fit1, corstr="AR-M-dependent")
fit3 <- update(fit1, corstr="Stationary-M-dependent(2)")
fit4 <- update(fit1, corstr="Exchangeable")
AGPC(fit1, fit2, fit3, fit4)

###### Example 2: Treatment for severe postnatal depression
data(depression)
mod2 <- depressd ~ visit + group
fit1 <- glmgee(mod2, id=subj, family=binomial(logit), data=depression)
fit2 <- update(fit1, corstr="AR-M-dependent")
fit3 <- update(fit1, corstr="Stationary-M-dependent(2)")
fit4 <- update(fit1, corstr="Exchangeable")
AGPC(fit1, fit2, fit3, fit4)

###### Example 3: Treatment for severe postnatal depression (2)
mod3 <- dep ~ visit*group
fit1 <- glmgee(mod3, id=subj, family=gaussian(identity), data=depression)
fit2 <- update(fit1, corstr="AR-M-dependent")
fit3 <- update(fit1, corstr="Exchangeable")
AGPC(fit1, fit2, fit3)

Description

A total of 1151 women completed menstrual diaries. The diary data were used to generate a binary sequence for each woman, indicating whether or not she had experienced amenorrhea (the absence of menstrual bleeding for a specified number of days) on the day of randomization and three additional 90-day intervals. The goal of this trial was to compare the two treatments (100 mg or 150 mg of depot-medroxyprogesterone acetate (DMPA)) in terms of how the rates of amenorrhea change over time with continued use of the contraceptive method.

Usage

data(amenorrhea)

Format

A data frame with 4604 rows and 4 variables:

ID

a numeric vector indicating the woman's ID.

Dose

a factor with two levels: "100mg" for treatment with 100 mg injection; and "150mg" for treatment with 150 mg injection.

Time

a numeric vector indicating the number of 90-day intervals since the trial beagn.

amenorrhea

a numeric vector indicating the amenorrhea status (1 for amenorrhea; 0 otherwise).

References

Fitzmaurice G.M., Laird N.M., Ware J.H. (2011, page 397). Applied Longitudinal Analysis. 2nd ed. John Wiley & Sons.

Machin D., Farley T.M., Busca B., Campbell M.J., d'Arcangues C. (1988) Assessing changes in vaginal bleeding patterns in contracepting women. Contraception, 38, 165-79.

Examples

data(amenorrhea)
dev.new()
amenorrhea2 <- aggregate(amenorrhea ~ Time + Dose,mean,data=amenorrhea,na.rm=TRUE)
barplot(100*amenorrhea ~ Dose+Time,data=amenorrhea2,beside=TRUE,col=c("blue","yellow"),ylab="%")
legend("topleft",legend=c("100 mg","150 mg"),fill=c("blue","yellow"),title="Dose",bty="n")

Description

Allows to compare nested generalized estimating equations using the Wald and generalized score tests.

Usage

## S3 method for class 'glmgee'
anova(
  object,
  ...,
  test = c("wald", "score"),
  verbose = TRUE,
  varest = c("robust", "df-adjusted", "model", "bias-corrected")
)

Arguments

Value

A matrix with three columns which contains the following:

Chi

The value of the statistic of the test.

df

The number of degrees of freedom.

Pr(>Chi)

The p-value of the test computed using the Chi-square distribution.

References

Rotnitzky A., Jewell P. (1990) Hypothesis Testing of Regression Parameters in Semiparametric Generalized Linear Models for Cluster Correlated Data. Biometrika 77:485-497.

Boos D.D. (1992) On Generalized Score Tests. The American Statistician 46:327-333.

Vanegas L.H., Rondon L.M., Paula G.A. (2023) Generalized Estimating Equations using the new R package glmtoolbox. The R Journal 15:105-133.

Boos D. (1992) On Generalized Score Tests. American Statistician 46:327–33.

Rotnitzky A., Jewell N.P. (1990). Hypothesis Testing of Regression Parameters in Semiparametric Generalized Linear Models for Cluster Correlated Data. Biometrika 77:485-497.

Examples

###### Example 1: Effect of ozone-enriched atmosphere on growth of sitka spruces
data(spruces)
mod <- size ~ poly(days,4)
fit1 <- glmgee(mod, id=tree, family=Gamma(log), data=spruces, corstr="AR-M-dependent")
fit2 <- update(fit1, . ~ . + treat)
fit3 <- update(fit2, . ~ . + poly(days,4):treat)
anova(fit1,fit2,fit3,test="wald")
anova(fit3,test="wald")

###### Example 2: Treatment for severe postnatal depression
data(depression)
mod2 <- depressd ~ group
fit1 <- glmgee(mod2, id=subj, family=binomial(logit), corstr="AR-M-dependent", data=depression)
fit2 <- update(fit1, . ~ . + visit)
fit3 <- update(fit2, . ~ . + group:visit)
anova(fit1,fit2,fit3,test="score")
anova(fit3,test="score")

Description

Allows to use the likelihood-ratio test to compare nested models in generalized nonlinear models.

Usage

## S3 method for class 'gnm'
anova(object, ..., verbose = TRUE)

Arguments

Value

A matrix with the following three columns:

Examples

###### Example: The effects of fertilizers on coastal Bermuda grass
data(Grass)
fit1 <- gnm(Yield ~ b0 + b1/(Nitrogen + a1) + b2/(Phosphorus + a2) + b3/(Potassium + a3),
            family=gaussian(inverse), start=c(b0=0.1,b1=13,b2=1,b3=1,a1=45,a2=15,a3=30), data=Grass)

fit2 <- update(fit1,Yield ~ I(b0 + b2/(Phosphorus + a2) + b3/(Potassium + a3)),
               start=c(b0=0.1,b2=1,b3=1,a2=15,a3=30))

anova(fit2,fit1)

Description

Allows to compare nested models for regression models based on the negative binomial, beta-binomial, and random-clumped binomial distributions, which are alternatives to the Poisson and binomial regression models under the presence of overdispersion. The comparisons are performed by using the Wald, score, gradient or likelihood ratio tests.

Usage

## S3 method for class 'overglm'
anova(object, ..., test = c("wald", "lr", "score", "gradient"), verbose = TRUE)

Arguments

Value

A matrix with the following three columns:

References

Buse A. (1982) The Likelihood Ratio, Wald, and Lagrange Multiplier Tests: An Expository Note. The American Statistician 36, 153-157.

Terrell G.R. (2002) The gradient statistic. Computing Science and Statistics 34, 206–215.

Examples

## Example 1: Self diagnozed ear infections in swimmers
data(swimmers)
fit1 <- overglm(infections ~ frequency + location + age + gender, family="nb1(log)", data=swimmers)
anova(fit1, test="wald")
anova(fit1, test="score")
anova(fit1, test="lr")
anova(fit1, test="gradient")

## Example 2: Agents to stimulate cellular differentiation
data(cellular)
fit2 <- overglm(cbind(cells,200-cells) ~ tnf*ifn, family="bb(logit)", data=cellular)
anova(fit2, test="wald")
anova(fit2, test="score")
anova(fit2, test="lr")
anova(fit2, test="gradient")

Description

Allows to compare nested models for regression models used to deal with zero-excess in count data. The comparisons are performed by using the Wald, score, gradient or likelihood ratio tests.

Usage

## S3 method for class 'zeroinflation'
anova(
  object,
  ...,
  test = c("wald", "lr", "score", "gradient"),
  verbose = TRUE,
  submodel = c("counts", "zeros")
)

Arguments

Value

A matrix with the following three columns:

Chi

The value of the statistic of the test,

Df

The number of degrees of freedom,

Pr(>Chi)

The p-value of the test test computed using the Chi-square distribution.

References

Buse A. (1982) The Likelihood Ratio, Wald, and Lagrange Multiplier Tests: An Expository Note. The American Statistician 36, 153-157.

Terrell G.R. (2002) The gradient statistic. Computing Science and Statistics 34, 206–215.

Examples

####### Example 1: Article production by graduate students in biochemistry PhD programs
bioChemists <- pscl::bioChemists
fit1 <- zeroinf(art ~ fem + kid5 + ment | ment, family="nb1(log)", data = bioChemists)
anova(fit1,test="wald")
anova(fit1,test="lr")
anova(fit1,test="score")
anova(fit1,test="gradient")

fit2 <- zeroalt(art ~ fem + kid5 + ment, family="nb1(log)", data = bioChemists)
anova(fit2,submodel="zeros",test="wald")
anova(fit2,submodel="zeros",test="lr")
anova(fit2,submodel="zeros",test="score")
anova(fit2,submodel="zeros",test="gradient")

Description

Allows to compare nested generalized linear models using Wald, score, gradient, and likelihood ratio tests.

Usage

anova2(
  object,
  ...,
  test = c("wald", "lr", "score", "gradient"),
  verbose = TRUE
)

Arguments

Details

The Wald, Rao's score and Terrell's gradient tests are performed using the expected Fisher information matrix.

Value

A matrix with three columns which contains the following:

Chi

The value of the statistic of the test.

Df

The number of degrees of freedom.

Pr(>Chi)

The p-value of the test computed using the Chi-square distribution.

References

Buse A. (1982) The Likelihood Ratio, Wald, and Lagrange Multiplier Tests: An Expository Note. The American Statistician 36, 153-157.

Terrell G.R. (2002) The gradient statistic. Computing Science and Statistics 34, 206 – 215.

Examples

## Example 1
Auto <- ISLR::Auto
fit1 <- glm(mpg ~ weight, family=inverse.gaussian("log"), data=Auto)
fit2 <- update(fit1, . ~ . + horsepower)
fit3 <- update(fit2, . ~ . + horsepower:weight)
anova2(fit1, fit2, fit3, test="lr")
anova2(fit1, fit2, fit3, test="score")
anova2(fit1, fit2, fit3, test="wald")
anova2(fit1, fit2, fit3, test="gradient")

## Example 2
burn1000 <- aplore3::burn1000
mod <- death ~ age + tbsa + inh_inj
fit1 <- glm(mod, family=binomial("logit"), data=burn1000)
fit2 <- update(fit1, . ~ . + inh_inj + age*inh_inj + tbsa*inh_inj)
anova2(fit1, fit2, test="lr")
anova2(fit1, fit2, test="score")
anova2(fit1, fit2, test="wald")
anova2(fit1, fit2, test="gradient")

## Example 3
data(aucuba)
fit <- glm(lesions ~ 1 + time, family=poisson("log"), data=aucuba)
anova2(fit, test="lr")
anova2(fit, test="score")
anova2(fit, test="wald")
anova2(fit, test="gradient")

Description

The investigators counted the number of lesions of Aucuba mosaic virus developing after exposure to X rays for various times.

Usage

data(aucuba)

Format

A data frame with 7 rows and 2 variables:

time

a numeric vector giving the minutes of exposure.

lesions

a numeric vector giving the counts of lesions, in hundreds.

References

Snedecor G.W., Cochran W.G. (1989) Statistical Methods, Eight Edition, Iowa State University Press, Ames.

Examples

data(aucuba)
dev.new()
barplot(lesions ~ time, col="red", data=aucuba)

Description

Best subset selection by exhaustive search in generalized linear models.

Usage

bestSubset(
  object,
  nvmax = 8,
  nbest = 1,
  force.in = NULL,
  force.out = NULL,
  verbose = TRUE,
  digits = max(3, getOption("digits") - 2)
)

Arguments

Details

In order to apply the "best subset" selection, an exhaustive search is conducted, separately for every size from $i$ to nvmax, to identify the model with the smallest deviance value. Therefore, if, for a fixed model size, the interest model selection criteria reduce to monotone functions of deviance, thus differing only in the way the sizes of the models are compared, then the results of the "best subset" selection do not depend upon the choice of the trade-off between goodness-of-fit and complexity on which they are based.

Examples

###### Example 1: Fuel consumption of automobiles
Auto <- ISLR::Auto
Auto2 <- within(Auto, origin <- factor(origin))
mod <- mpg ~ cylinders + displacement + acceleration + origin + horsepower*weight
fit1 <- glm(mod, family=inverse.gaussian(log), data=Auto2)
out1 <- bestSubset(fit1)
out1

###### Example 2: Patients with burn injuries
burn1000 <- aplore3::burn1000
burn1000 <- within(burn1000, death <- factor(death, levels=c("Dead","Alive")))
mod <- death ~ gender + race + flame + age*tbsa*inh_inj
fit2 <- glm(mod, family=binomial(logit), data=burn1000)
out2 <- bestSubset(fit2)
out2

###### Example 3: Advertising
data(advertising)
fit3 <- glm(sales ~ log(TV)*radio*newspaper, family=gaussian(log), data=advertising)
out3 <- bestSubset(fit3)
out3

Description

Female mice were continuously fed dietary concentrations of 2-Acetylaminofluorene (2-AAF), a carcinogenic and mutagenic derivative of fluorene. Serially sacrificed, dead or moribund mice were examined for tumors and deaths dates recorded. These data consist of the incidences of bladder neoplasms in mice observed during 33 months.

Usage

data(bladder)

Format

A data frame with 8 rows and 3 variables:

dose

a numeric vector giving the dose, in parts per $10^4$, of 2-AAF.

exposed

a numeric vector giving the number of mice exposed to each dose of 2-AAF.

cancer

a numeric vector giving the number of mice with bladder cancer for each dose of 2-AAF.

References

Zhang H., Zelterman D. (1999) Binary Regression for Risks in Excess of Subject-Specific Thresholds. Biometrics 55:1247-1251.

See Also

liver

Examples

data(bladder)
dev.new()
barplot(100*cancer/exposed ~ dose, beside=TRUE, data=bladder, col="red",
        xlab="Dose of 2-AAF", ylab="% of mice with bladder cancer")

Description

Computes the Box-Tidwell power transformations of the predictors in a regression model.

Usage

BoxTidwell(
  object,
  transf,
  epsilon = 1e-04,
  maxiter = 30,
  trace = FALSE,
  digits = max(3, getOption("digits") - 2),
  ...
)

Arguments

Value

Two matrices with the values of marginal and omnibus tests.

Description

Computes the Box-Tidwell power transformations of the predictors in a generalized linear model.

Usage

## S3 method for class 'glm'
BoxTidwell(
  object,
  transf,
  epsilon = 1e-04,
  maxiter = 30,
  trace = FALSE,
  digits = max(3, getOption("digits") - 2),
  ...
)

Arguments

Value

a list list with components including

References

Box G.E.P., Tidwell P.W. (1962) Transformation of the independent variables. Technometrics 4, 531-550.

Fox J. (2016) Applied Regression Analysis and Generalized Linear Models, Third Edition. Sage.

See Also

BoxTidwell.lm

Examples

###### Example 1: Skin cancer in women
data(skincancer)
fit1 <- glm(cases ~ age + city, offset=log(population), family=poisson(log), data=skincancer)
AIC(fit1)
BoxTidwell(fit1, transf= ~ age)
fit1 <- update(fit1,formula=. ~ I(age^(-1/2)) + city)
AIC(fit1)

###### Example 3: Gas mileage
data(Auto, package="ISLR")
fit3 <- glm(mpg ~ horsepower + weight, family=inverse.gaussian(log), data=Auto)
AIC(fit3)
BoxTidwell(fit3, transf= ~ horsepower + weight)
fit3 <- update(fit3,formula=. ~ I(horsepower^(-1/3)) + weight)
AIC(fit3)

###### Example 4: Advertising
data(advertising)
fit4 <- glm(sales ~ TV + radio, family=gaussian(log), data=advertising)
AIC(fit4)
BoxTidwell(fit4, transf= ~ TV)
fit4 <- update(fit4,formula=. ~ I(TV^(1/10)) + radio)
AIC(fit4)

###### Example 5: Leukaemia Patients
data(leuk, package="MASS")
fit5 <- glm(ifelse(time>=52,1,0) ~ ag + wbc, family=binomial, data=leuk)
AIC(fit5)
BoxTidwell(fit5, transf= ~ wbc)
fit5 <- update(fit5,formula=. ~ ag + I(wbc^(-0.18)))
AIC(fit5)

Description

Computes the Box-Tidwell power transformations of the predictors in a normal linear model.

Usage

## S3 method for class 'lm'
BoxTidwell(
  object,
  transf,
  epsilon = 1e-04,
  maxiter = 30,
  trace = FALSE,
  digits = max(3, getOption("digits") - 2),
  ...
)

Arguments

Value

a list list with components including

References

Box G.E.P., Tidwell P.W. (1962) Transformation of the independent variables. Technometrics 4, 531-550.

Fox J. (2016) Applied Regression Analysis and Generalized Linear Models, Third Edition. Sage.

See Also

BoxTidwell.glm

Examples

###### Example 1: Hill races in Scotland
data(races)
fit1 <- lm(rtime ~ distance + cclimb, data=races)
AIC(fit1)
BoxTidwell(fit1, transf= ~ distance + cclimb)
fit1 <- update(fit1,formula=rtime ~ distance + I(cclimb^2))
AIC(fit1)

###### Example 2: Gasoline yield
fit2 <- lm(mpg ~ hp + wt + am, data=mtcars)
AIC(fit2)
BoxTidwell(fit2, transf= ~ hp + wt)
fit2 <- update(fit2,formula=mpg ~ log(hp) + log(wt) + am)
AIC(fit2)

###### Example 3: New York Air Quality Measurements
fit3 <- lm(log(Ozone) ~ Solar.R + Wind + Temp, data=airquality)
AIC(fit3)
BoxTidwell(fit3, transf= ~ Solar.R + Wind + Temp)
fit3 <- update(fit3,formula=log(Ozone) ~ log(Solar.R) + Wind + Temp)
AIC(fit3)

###### Example 4: Heat capacity of hydrobromic acid
data(heatcap,package="GLMsData")
fit4 <- lm(Cp ~ Temp, data=heatcap)
AIC(fit4)
BoxTidwell(fit4, transf= ~ Temp)
fit4 <- update(fit4,formula=Cp ~ I(Temp^5))
AIC(fit4)

###### Example 5: Age and Eye Lens Weight of Rabbits in Australia
data(rabbits)
fit5 <- lm(log(wlens) ~ age, data=rabbits)
AIC(fit5)
BoxTidwell(fit5, transf= ~ age)
fit5 <- update(fit5,formula=log(wlens) ~ I(age^(-1/3)))
AIC(fit5)

Description

These data correspond to the (average) body weight and the (average) brain weight for sixty-two species of mammals.

Usage

data(brains)

Format

A data frame with 62 rows and 3 variables:

Specie

a character string giving the species name.

BrainWt

a numeric vector indicating the average brain weight, in grams.

BodyWt

a numeric vector indicating the average body weight, in kilograms.

References

Allison T., Cicchetti D. (1976). Sleep in mammals: Ecology and constitutional correlates. Science 194:732-734.

Weisberg S. (2005). Applied Linear Regression, 3rd edition. Wiley, New York.

Examples

data(brains)
xlab <- "log(Body Weight)"
ylab <- "log(Brain Weight)"
dev.new()
with(brains,plot(log(BodyWt),log(BrainWt),pch=20,xlab=xlab,ylab=ylab))

Description

Data on technical support calls after a product release. Using this information, new products could be allocated technical support resources.

Usage

data(calls)

Format

A data frame with 16 rows and 2 variables:

week

a numeric vector indicating number of weeks that have elapsed since the product’s release.

calls

a numeric vector indicating the number of technical support calls.

Source

https://documentation.sas.com/doc/en/pgmsascdc/9.4_3.3/statug/statug_mcmc_examples12.htm

Examples

data(calls)
dev.new()
with(calls,plot(week,calls,xlab="The number of weeks since the release of the product",
     pch=16,col="blue",ylab="Technical support calls"))

Description

In a biomedical study of the immuno-activating ability of two agents, TNF (tumor necrosis factor) and IFN (interferon), to induce cell differentiation, the number of cells that exhibited markers of differentiation after exposure to TNF and IFN was recorded. At each of the 16 TNF/INF dose combinations, 200 cells were examined. The main question is whether the two agents stimulate cell differentiation synergistically or independently.

Usage

data(cellular)

Format

A data frame with 16 rows and 3 variables:

cells

a numeric vector giving the number of cells that exhibited markers of differentiation after exposure to the dose of the two agents

tnf

a numeric vector giving the dose (U/ml) of TNF

ifn

a numeric vector giving the dose (U/ml) of IFN

References

Piegorsch W.W., Weinberg C.R., Margolin B.H. (1988) Exploring simple independent action in multifactor tables of proportions. Biometrics 44:595-603.

Vanegas L.H., Rondon L.M. (2020) A data transformation to deal with constant under/over-dispersion in binomial and poisson regression models. Journal of Statistical Computation and Simulation 90:1811-1833.

Examples

data(cellular)
dev.new()
barplot(100*cells/200 ~ ifn + tnf, beside=TRUE, data=cellular, col=terrain.colors(4),
        xlab="Dose of TNF", ylab="% of cells with markers of differentiation")
legend("topleft", legend=c("0","4","20","100"), fill=terrain.colors(4),
       title="Dose of IFN", bty="n")

Description

Inflation of the abdomen during laparoscopic cholecystectomy (removal of the gallbladder) separates the liver from the diaphragm and strains the attachments that connect both. This strain is felt as a referred shoulder pain. Suction to remove residual gas may reduce shoulder pain. There were 22 subjects randomized in the active group (with abdominal suction) and 19 subjects randomized in the control group (without abdominal suction). After laparoscopic surgery, patients were asked to rate their shoulder pain on a visual analog scale morning and afternoon for three days after the operation (a total of six different times). The scale was coded into five ordered categories where a pain score of 1 indicated "low pain" and a score of 5 reflected "high pain".

Usage

data(cholecystectomy)

Format

A data frame with 246 rows and 7 variables:

id

a numeric vector with the identifier of the patient.

treatment

a factor indicating the treatment received by the patient: abdominal suction ("A") and placebo ("P").

gender

a factor indicating the gender of the patient: female ("F") and male ("M").

age

a numeric vector indicating the age of the patient, in years.

time

a numeric vector indicating the occasion the patient was asked to rate their shoulder pain after the laparoscopic surgery: integers from 1 to 6.

pain

a numeric vector indicating the shoulder pain rated by the patient on a scale coded into five ordered categories, where 1 indicated "low pain" and 5 reflected "high pain".

pain2

a numeric vector indicating the shoulder pain rated by the patient and coded as 1 for the two first categories of pain and 0 for other cases.

References

Jorgensen J.O., Gillies R.B., Hunt D.R., Caplehorn J.R.M., Lumley T. (1995) A simple and effective way to reduce postoperative pain after laparoscopic cholecystectomy. Australian and New Zealand Journal of Surgery 65:466–469.

Lumley T. (1996) Generalized Estimating Equations for Ordinal Data: A Note on Working Correlation Structures. Biometrics 52:354–361.

Morel J.G., Nagaraj N.K. (2012) Overdispersion Models in SAS. SAS Institute Inc., Cary, North Carolina, USA.

Examples

data(cholecystectomy)
out <- aggregate(pain2 ~ treatment + time, data=cholecystectomy, mean)
dev.new()
barplot(100*pain2 ~ treatment + time, beside=TRUE, data=out, xlab="Time",
        col=c("yellow","blue"), ylab="% of patients with low pain")
legend("topleft", c("Placebo","Abdominal suction"), fill=c("yellow","blue"),
       title="Treatment", cex=0.9, bty="n")

Description

Computes the Correlation Information Criterion (CIC) for one or more objects of the class glmgee.

Usage

CIC(..., verbose = TRUE, digits = max(3, getOption("digits") - 2))

Arguments

Value

A data.frame with the values of the CIC for each glmgee object in the input.

References

Hin L.-Y., Wang Y.-G. (2009) Working-Correlation-Structure Identification in Generalized Estimating Equations. Statistics in Medicine, 28:642-658.

Hin L.-Y., Carey V.J., Wang Y.-G. (2007) Criteria for Working–Correlation–Structure Selection in GEE: Assessment via Simulation. The American Statistician 61:360–364.

Vanegas L.H., Rondon L.M., Paula G.A. (2023) Generalized Estimating Equations using the new R package glmtoolbox. The R Journal 15:105-133.

See Also

QIC, GHYC, RJC, AGPC, SGPC

Examples

###### Example 1: Effect of ozone-enriched atmosphere on growth of sitka spruces
data(spruces)
mod1 <- size ~ poly(days,4) + treat
fit1 <- glmgee(mod1, id=tree, family=Gamma(log), data=spruces)
fit2 <- update(fit1, corstr="AR-M-dependent")
fit3 <- update(fit1, corstr="Stationary-M-dependent(2)")
fit4 <- update(fit1, corstr="Exchangeable")
CIC(fit1, fit2, fit3, fit4)

###### Example 2: Treatment for severe postnatal depression
data(depression)
mod2 <- depressd ~ visit + group
fit1 <- glmgee(mod2, id=subj, family=binomial(logit), data=depression)
fit2 <- update(fit1, corstr="AR-M-dependent")
fit3 <- update(fit1, corstr="Stationary-M-dependent(2)")
fit4 <- update(fit1, corstr="Exchangeable")
CIC(fit1, fit2, fit3, fit4)

###### Example 3: Treatment for severe postnatal depression (2)
mod3 <- dep ~ visit*group
fit1 <- glmgee(mod3, id=subj, family=gaussian(identity), data=depression)
fit2 <- update(fit1, corstr="AR-M-dependent")
fit3 <- update(fit1, corstr="Exchangeable")
CIC(fit1, fit2, fit3)

Description

Computes confidence intervals based on Wald test for a generalized nonlinear model.

Usage

## S3 method for class 'gnm'
confint(
  object,
  parm,
  level = 0.95,
  contrast,
  digits = max(3, getOption("digits") - 2),
  dispersion = NULL,
  verbose = TRUE,
  ...
)

Arguments

Details

The approximate 100(level)% confidence interval for $\beta$ based on the Wald test.

Value

A matrix with so many rows as parameters in the "linear" predictor and two columns: "Lower limit" and "Upper limit".

Examples

###### Example 1: The effects of fertilizers on coastal Bermuda grass
data(Grass)
fit1 <- gnm(Yield ~ b0 + b1/(Nitrogen + a1) + b2/(Phosphorus + a2) + b3/(Potassium + a3),
            family=gaussian(inverse), start=c(b0=0.1,b1=13,b2=1,b3=1,a1=45,a2=15,a3=30), data=Grass)

confint(fit1, level=0.95)

###### Example 2: Assay of an Insecticide with a Synergist
data(Melanopus)
fit2 <- gnm(Killed/Exposed ~ b0 + b1*log(Insecticide-a1) + b2*Synergist/(a2 + Synergist),
            family=binomial(logit), weights=Exposed, start=c(b0=-3,b1=1.2,a1=1.7,b2=1.7,a2=2),
		   data=Melanopus)

confint(fit2, level=0.95)

###### Example 3: Developmental rate of Drosophila melanogaster
data(Drosophila)
fit3 <- gnm(Duration ~ b0 + b1*Temp + b2/(Temp-a), family=Gamma(log),
            start=c(b0=3,b1=-0.25,b2=-210,a=55), weights=Size, data=Drosophila)

confint(fit3, level=0.95)

###### Example 4: Radioimmunological Assay of Cortisol
data(Cortisol)
fit4 <- gnm(Y ~ b0 + (b1-b0)/(1 + exp(b2+ b3*lDose))^b4, family=Gamma(identity),
            start=c(b0=130,b1=2800,b2=3,b3=3,b4=0.5), data=Cortisol)

### Confidence Interval for 'b1-b0'
confint(fit4, level=0.95, contrast=matrix(c(-1,1,0,0,0),1,5))

### Confidence Intervals for 'b0', 'b1', 'b2', 'b3', 'b4'
confint(fit4, level=0.95, contrast=diag(5))

Description

Computes confidence intervals based on Wald, likelihood-ratio, Rao's score or Terrell's gradient tests for a generalized linear model.

Usage

confint2(
  model,
  level = 0.95,
  test = c("wald", "lr", "score", "gradient"),
  digits = max(3, getOption("digits") - 2),
  verbose = TRUE
)

Arguments

Details

The approximate 100(level)% confidence interval for $\beta$ based on the test test is the set of values of $\beta_0$ for which the hypothesis $H_0$: $\beta=\beta_0$ versus $H_1$: $\beta!=\beta_0$ is not rejected at the approximate significance level of 100(1-level)%. The Wald, Rao's score and Terrell's gradient tests are performed using the expected Fisher information matrix.

Value

A matrix with so many rows as parameters in the linear predictor and two columns: "Lower limit" and "Upper limit".

References

Buse A. (1982) The Likelihood Ratio, Wald, and Lagrange Multiplier Tests: An Expository Note. The American Statistician 36, 153-157.

Terrell G.R. (2002) The gradient statistic. Computing Science and Statistics 34, 206 – 215.

Examples

###### Example 1: Fuel consumption of automobiles
Auto <- ISLR::Auto
fit1 <- glm(mpg ~ weight*horsepower, family=inverse.gaussian("log"), data=Auto)
confint2(fit1, test="lr")
confint2(fit1, test="score")

###### Example 2: Patients with burn injuries
burn1000 <- aplore3::burn1000
burn1000 <- within(burn1000, death <- factor(death, levels=c("Dead","Alive")))
fit2 <- glm(death ~ age*inh_inj + tbsa*inh_inj, family=binomial("logit"), data=burn1000)
confint2(fit2, test="lr")
confint2(fit2, test="gradient")

Description

Produces an approximation, better known as the one-step aproximation, of the Cook's distance, which is aimed to measure the effect on the estimates of the parameters in the linear predictor of deleting each cluster/observation in turn. This function also can produce a cluster/observation-index plot of the Cook's distance for all parameters in the linear predictor or for some subset of them (via the argument coefs).

Usage

## S3 method for class 'glmgee'
cooks.distance(
  model,
  method = c("Preisser-Qaqish", "full"),
  level = c("clusters", "observations"),
  plot.it = FALSE,
  coefs,
  identify,
  varest = c("robust", "df-adjusted", "model", "bias-corrected"),
  ...
)

Arguments

Details

The Cook's distance consists of the distance between two estimates of the parameters in the linear predictor using a metric based on the (estimate of the) variance-covariance matrix. For the cluster-level, the first one set of estimates is computed from a dataset including all clusters/observations, and the second one is computed from a dataset in which the i-th cluster is excluded. To avoid computational burden, the second set of estimates is replaced by its one-step approximation. See the dfbeta.glmgee documentation.

Value

A matrix as many rows as clusters/observations in the sample and one column with the values of the Cook's distance.

References

Pregibon D. (1981). Logistic regression diagnostics. The Annals of Statistics 9, 705-724.

Preisser J.S., Qaqish B.F. (1996) Deletion diagnostics for generalised estimating equations. Biometrika 83:551–562.

Hammill B.G., Preisser J.S. (2006) A SAS/IML software program for GEE and regression diagnostics. Computational Statistics & Data Analysis 51:1197-1212.

Vanegas L.H., Rondon L.M., Paula G.A. (2023) Generalized Estimating Equations using the new R package glmtoolbox. The R Journal 15:105-133.

Examples

###### Example 1: Effect of ozone-enriched atmosphere on growth of sitka spruces
data(spruces)
mod1 <- size ~ poly(days,4) + treat
fit1 <- glmgee(mod1, id=tree, family=Gamma(log), data=spruces, corstr="AR-M-dependent")

### Cook's distance for all parameters in the linear predictor
cooks.distance(fit1, method="full", plot.it=TRUE, col="red", lty=1, lwd=1, cex=0.8,
               col.lab="blue", col.axis="blue", col.main="black", family="mono")

### Cook's distance for the parameter associated to the variable 'treat'
cooks.distance(fit1, coef="treat", method="full", plot.it=TRUE, col="red", lty=1,
               lwd=1, col.lab="blue", col.axis="blue", col.main="black", cex=0.8)

###### Example 2: Treatment for severe postnatal depression
data(depression)
mod2 <- depressd ~ visit + group
fit2 <- glmgee(mod2, id=subj, family=binomial(logit), corstr="AR-M-dependent", data=depression)

### Cook's distance for all parameters in the linear predictor
cooks.distance(fit2, method="full", plot.it=TRUE, col="red", lty=1, lwd=1, cex=0.8,
               col.lab="blue", col.axis="blue", col.main="black", family="mono")

### Cook's distance for the parameter associated to the variable 'group'
cooks.distance(fit2, coef="group", method="full", plot.it=TRUE, col="red", lty=1,
               lwd=1, col.lab="blue", col.axis="blue", col.main="black", cex=0.8)

Description

Produces an approximation of the Cook's distance, better known as the one-step approximation, for measuring the effect of deleting each observation in turn on the estimates of the parameters in a linear predictor. Additionally, this function can produce an index plot of Cook's distance for all or a subset of the parameters in the linear predictor (via the argument coefs).

Usage

## S3 method for class 'gnm'
cooks.distance(model, plot.it = FALSE, dispersion = NULL, coefs, identify, ...)

Arguments

Details

The Cook's distance consists of the distance between two estimates of the parameters in the linear predictor using a metric based on the (estimate of the) variance-covariance matrix. The first one set of estimates is computed from a dataset including all individuals, and the second one is computed from a dataset in which the i-th individual is excluded. To avoid computational burden, the second set of estimates is replaced by its one-step approximation. See the dfbeta.overglm documentation.

Value

A matrix as many rows as individuals in the sample and one column with the values of the Cook's distance.

Examples

###### Example 1: The effects of fertilizers on coastal Bermuda grass
data(Grass)
fit1 <- gnm(Yield ~ b0 + b1/(Nitrogen + a1) + b2/(Phosphorus + a2) + b3/(Potassium + a3),
            family=gaussian(inverse), start=c(b0=0.1,b1=13,b2=1,b3=1,a1=45,a2=15,a3=30), data=Grass)

cooks.distance(fit1, plot.it=TRUE, col="red", lty=1, lwd=1,
  col.lab="blue", col.axis="blue", col.main="black", family="mono", cex=0.8)

###### Example 2: Assay of an Insecticide with a Synergist
data(Melanopus)
fit2 <- gnm(Killed/Exposed ~ b0 + b1*log(Insecticide-a1) + b2*Synergist/(a2 + Synergist),
            family=binomial(logit), weights=Exposed, start=c(b0=-3,b1=1.2,a1=1.7,b2=1.7,a2=2),
		   data=Melanopus)

### Cook's distance just for the parameter "b1"
cooks.distance(fit2, plot.it=TRUE, coef="b1", col="red", lty=1, lwd=1,
  col.lab="blue", col.axis="blue", col.main="black", family="mono", cex=0.8)

###### Example 3: Developmental rate of Drosophila melanogaster
data(Drosophila)
fit3 <- gnm(Duration ~ b0 + b1*Temp + b2/(Temp-a), family=Gamma(log),
            start=c(b0=3,b1=-0.25,b2=-210,a=55), weights=Size, data=Drosophila)

cooks.distance(fit3, plot.it=TRUE, col="red", lty=1, lwd=1,
  col.lab="blue", col.axis="blue", col.main="black", family="mono", cex=0.8)

###### Example 4: Radioimmunological Assay of Cortisol
data(Cortisol)
fit4 <- gnm(Y ~ b0 + (b1-b0)/(1 + exp(b2+ b3*lDose))^b4, family=Gamma(identity),
            start=c(b0=130,b1=2800,b2=3,b3=3,b4=0.5), data=Cortisol)

cooks.distance(fit4, plot.it=TRUE, col="red", lty=1, lwd=1,
  col.lab="blue", col.axis="blue", col.main="black", family="mono", cex=0.8)

Description

Produces an approximation, better known as the one-step approximation, of the Cook's distance, which is aimed to measure the effect on the estimates of the parameters in the linear predictor of deleting each observation in turn. This function also can produce an index plot of the Cook's distance for all parameters in the linear predictor or for some subset of them (via the argument coefs).

Usage

## S3 method for class 'overglm'
cooks.distance(model, plot.it = FALSE, coefs, identify, ...)

Arguments

Details

The Cook's distance consists of the distance between two estimates of the parameters in the linear predictor using a metric based on the (estimate of the) variance-covariance matrix. The first one set of estimates is computed from a dataset including all individuals, and the second one is computed from a dataset in which the i-th individual is excluded. To avoid computational burden, the second set of estimates is replaced by its one-step approximation. See the dfbeta.overglm documentation.

Value

A matrix as many rows as individuals in the sample and one column with the values of the Cook's distance.

Examples

###### Example 1: Self diagnozed ear infections in swimmers
data(swimmers)
fit1 <- overglm(infections ~ frequency + location, family="nb1(log)", data=swimmers)

### Cook's distance for all parameters in the linear predictor
cooks.distance(fit1, plot.it=TRUE, col="red", lty=1, lwd=1, col.lab="blue",
               col.axis="blue", col.main="black", family="mono", cex=0.8)

### Cook's distance just for the parameter associated with 'frequency'
cooks.distance(fit1, plot.it=TRUE, coef="frequency", col="red", lty=1, lwd=1,
   col.lab="blue", col.axis="blue", col.main="black", family="mono", cex=0.8)

###### Example 2: Article production by graduate students in biochemistry PhD programs
bioChemists <- pscl::bioChemists
fit2 <- overglm(art ~ fem + kid5 + ment, family="nb1(log)", data = bioChemists)

### Cook's distance for all parameters in the linear predictor
cooks.distance(fit2, plot.it=TRUE, col="red", lty=1, lwd=1, col.lab="blue",
               col.axis="blue", col.main="black", family="mono", cex=0.8)

### Cook's distance just for the parameter associated with 'fem'
cooks.distance(fit2, plot.it=TRUE, coef="fem", col="red", lty=1, lwd=1,
   col.lab="blue", col.axis="blue", col.main="black", family="mono", cex=0.8)

###### Example 3: Agents to stimulate cellular differentiation
data(cellular)
fit3 <- overglm(cbind(cells,200-cells) ~ tnf + ifn, family="bb(logit)", data=cellular)

### Cook's distance for all parameters in the linear predictor
cooks.distance(fit3, plot.it=TRUE, col="red", lty=1, lwd=1, col.lab="blue",
               col.axis="blue", col.main="black", family="mono", cex=0.8)

### Cook's distance just for the parameter associated with 'tnf'
cooks.distance(fit3, plot.it=TRUE, coef="tnf", col="red", lty=1, lwd=1,
  col.lab="blue", col.axis="blue", col.main="black", family="mono", cex=0.8)

Description

Produces an approximation, better known as the one-step approximation, of the Cook's distance, which is aimed to measure the effect on the estimates of the parameters in the linear predictor of deleting each observation in turn. This function also can produce an index plot of the Cook's distance for all parameters in the linear predictor or for some subset of them (via the argument coefs).

Usage

## S3 method for class 'zeroinflation'
cooks.distance(
  model,
  submodel = c("counts", "zeros", "full"),
  plot.it = FALSE,
  coefs,
  identify,
  ...
)

Arguments

Details

The Cook's distance consists of the distance between two estimates of the parameters in the linear predictor using a metric based on the (estimate of the) variance-covariance matrix. The first one set of estimates is computed from a dataset including all individuals, and the second one is computed from a dataset in which the i-th individual is excluded. To avoid computational burden, the second set of estimates is replaced by its one-step approximation. See the dfbeta.zeroinflation documentation.

Value

A matrix as many rows as individuals in the sample and one column with the values of the Cook's distance.

Examples

####### Example 1: Self diagnozed ear infections in swimmers
data(swimmers)
fit <- zeroinf(infections ~ frequency + location, family="nb1(log)", data=swimmers)

### Cook's distance for all parameters in the "counts" model
cooks.distance(fit, submodel="counts", plot.it=TRUE, col="red", lty=1, lwd=1,
         col.lab="blue", col.axis="blue", col.main="black", family="mono", cex=0.8)

### Cook's distance for all parameters in the "zeros" model
cooks.distance(fit, submodel="zeros", plot.it=TRUE, col="red", lty=1, lwd=1,
         col.lab="blue", col.axis="blue", col.main="black", family="mono", cex=0.8)

Description

The amount of hormone contained in a preparation cannot be measured directly, so estimating an unknown dose of hormone involves a two-step process. A calibration curve must first be established, then the curve must be inverted to determine the hormone dose. The calibration curve is estimated using a radioimmunological assay.

Usage

data(Cortisol)

Format

A data frame with 64 rows and 2 variables:

lDose

a numeric vector indicating the logarithm in base 10 of the dose.

Y

a numeric vector indicating the response, in counts per minute.

References

Huet S., Bouvier A., Poursat M.-A., Jolivet E. (2004). Statistical tools for nonlinear regression : a practical guide with S-PLUS and R examples. 2nd Edition. Springer, New York.

Examples

data(Cortisol)
dev.new()
with(Cortisol, plot(lDose,Y,xlab="Log10(Dose, in ng/0.1 ml)",pch=16,col="blue",
                    ylab="Response, in counts per minute"))

Description

The market research department of a soft drink manufacturer is investigating the effectiveness of a price discount coupon on the purchase of a two-litre beverage product. A sample of 5500 costumers received coupons for varying price discounts between 5 and 25 cents. The main objective of the analysis is to determine if the price discount affects the proportion of redeemed coupons after one month.

Usage

data(coupons)

Format

A data frame with 11 rows and 3 variables:

discounts

a numeric vector indicating the price discount, in cents.

costumers

a numeric vector indicating the number of customers who received coupons.

redeemed

a numeric vector indicating the number of redeemed coupons.

References

Montgomery D.C., Peck E.A., Vining G. (2012, page 464) Introduction to linear regression analysis. 5th ed. Berlin, Wiley.

Examples

dev.new()
data(coupons)
barplot(100*redeemed/costumers ~ discounts, data=coupons, xlab="Discount price",
        ylab="(%) Redeemed coupons", col="blue")

Description

These data arose from a study on the efficacy of oestrogen given transdermally for the treatment of severe postnatal depression. Women with major depression were randomly assigned to a placebo control group or an oestrogen patch group. Prior to the treatment all women were assessed by self-rated depressive symptoms on the Edinburgh Postnatal Depression Scale (EPDS). EPDS data were collected monthly for six months once treatment began. Higher EDPS scores are indicative of higher depression levels.

Usage

data(depression)

Format

A data frame with 427 rows and 5 variables:

subj

a numeric vector giving the identifier of each woman.

group

a factor giving the received treatment: "placebo" or "oestrogen".

visit

a numeric vector giving the number of months since the treatment began, where -1 indicates the pretreatment assessment of the EDPS.

dep

a numeric vector giving the value of the EDPS.

depressd

a numeric vector coded as 1 when the value of the EDPS is greater than or equal to 11 and coded as 0 in other cases.

Source

https://stats.oarc.ucla.edu/spss/library/spss-librarypanel-data-analysis-using-gee/

References

Gregoire A.J.P., Kumar R., Everitt B., Henderson A.F., Studd J.W.W. (1996) Transdermal oestrogen for treatment of severe postnatal depression, The Lancet 347:930-933.

Examples

data(depression)
dev.new()
boxplot(dep ~ visit, data=subset(depression,group=="placebo"), at=c(0:6) - 0.2,
        col="yellow", boxwex=0.3, xaxt="n", ylim=range(na.omit(depression$dep)),
        xlab="Months since the treatment began", ylab="EDPS")
boxplot(dep ~ visit, data=subset(depression,group=="oestrogen"), add=TRUE,
        at=c(0:6) + 0.2, col="blue", boxwex=0.3, xaxt="n")
axis(1, at=c(0:6), labels=c(-1,1:6))
legend("bottomleft", legend=c("placebo","oestrogen"), fill=c("yellow","blue"),
       title="Treatment", bty="n")

Description

Produces an approximation, better known as the one-step approximation, of the effect on the parameter estimates of deleting each cluster/observation in turn. This function also can produce an index plot of the Dfbeta Statistic for some parameters via the argument coefs.

Usage

## S3 method for class 'glmgee'
dfbeta(
  model,
  level = c("clusters", "observations"),
  method = c("Preisser-Qaqish", "full"),
  coefs,
  identify,
  ...
)

Arguments

Details

The one-step approximation (with the method "full") of the estimates of the parameters in the linear predictor of a GEE when the i-th cluster is excluded from the dataset is given by the vector obtained as the result of the first iteration of the fitting algorithm of that GEE when it is performed using: (1) a dataset in which the i-th cluster is excluded; and (2) a starting value which is the solution to the same GEE but based on the dataset inluding all clusters.

Value

A matrix with so many rows as clusters/observations in the sample and so many columns as parameters in the linear predictor. For clusters, the $i$-th row of that matrix corresponds to the difference between the estimates of the parameters in the linear predictor using all clustersand the one-step approximation of those estimates when the i-th cluster is excluded from the dataset.

References

Pregibon D. (1981). Logistic regression diagnostics. The Annals of Statistics 9, 705-724.

Preisser J.S., Qaqish B.F. (1996) Deletion diagnostics for generalised estimating equations. Biometrika 83:551–562.

Hammill B.G., Preisser J.S. (2006) A SAS/IML software program for GEE and regression diagnostics. Computational Statistics & Data Analysis 51:1197-1212.

Vanegas L.H., Rondon L.M., Paula G.A. (2023) Generalized Estimating Equations using the new R package glmtoolbox. The R Journal 15:105-133.

Examples

###### Example 1: Effect of ozone-enriched atmosphere on growth of sitka spruces
data(spruces)
mod1 <- size ~ poly(days,4) + treat
fit1 <- glmgee(mod1, id=tree, family=Gamma(log), corstr="AR-M-dependent", data=spruces)
dfbs1 <- dfbeta(fit1, method="full", coefs="treat", col="red", lty=1, lwd=1, col.lab="blue",
         col.axis="blue", col.main="black", family="mono", cex=0.8, main="treat")

### Calculation by hand of dfbeta for the tree labeled by "N1T01"
onestep1 <- glmgee(mod1, id=tree, family=Gamma(log), corstr="AR-M-dependent",
            data=spruces, start=coef(fit1), subset=c(tree!="N1T01"), maxit=1)

coef(fit1)-coef(onestep1)
dfbs1[rownames(dfbs1)=="N1T01",]

###### Example 2: Treatment for severe postnatal depression
data(depression)
mod2 <- depressd ~ visit + group
fit2 <- glmgee(mod2, id=subj, family=binomial(logit), corstr="AR-M-dependent",
               data=depression)

dfbs2 <- dfbeta(fit2, method="full", coefs="group" ,col="red", lty=1, lwd=1, col.lab="blue",
         col.axis="blue", col.main="black", family="mono", cex=0.8, main="group")

### Calculation by hand of dfbeta for the woman labeled by "18"
onestep2 <- glmgee(mod2, id=subj, family=binomial(logit), corstr="AR-M-dependent",
            data=depression, start=coef(fit2), subset=c(subj!=18), maxit=1)

coef(fit2)-coef(onestep2)
dfbs2[rownames(dfbs2)==18,]

###### Example 3: Treatment for severe postnatal depression (2)
mod3 <- dep ~ visit*group
fit3 <- glmgee(mod3, id=subj, family=gaussian(identity), corstr="AR-M-dependent",
               data=depression)

dfbs3 <- dfbeta(fit3, method="full", coefs="visit:group" ,col="red", lty=1, lwd=1, col.lab="blue",
         col.axis="blue", col.main="black", family="mono", cex=0.8, main="visit:group")

### Calculation by hand of dfbeta for the woman labeled by "18"
onestep3 <- glmgee(mod3, id=subj, family=gaussian(identity), corstr="AR-M-dependent",
            data=depression, start=coef(fit3), subset=c(subj!=18), maxit=1)

coef(fit3)-coef(onestep3)
dfbs3[rownames(dfbs3)==18,]

Description

Calculates an approximation of the parameter estimates that would be produced by deleting each case in turn, which is known as the one-step approximation. Additionally, the function can produce an index plot of the Dfbeta statistic for some parameter specified by the argument coefs.

Usage

## S3 method for class 'gnm'
dfbeta(model, coefs, identify, ...)

Arguments

Details

The one-step approximation of the parameters estimates when the $i$-th case is excluded from the dataset consists of the vector obtained as a result of the first iteration of the Fisher Scoring algorithm when it is performed using: (1) a dataset in which the $i$-th case is excluded; and (2) a starting value that is the estimate of the same model but based on the dataset including all cases.

Value

A matrix with as many rows as cases in the sample and as many columns as parameters in the linear predictor. The $i$-th row in that matrix corresponds to the difference between the parameters estimates obtained using all cases and the one-step approximation of those estimates when excluding the $i$-th case from the dataset.

References

Pregibon D. (1981). Logistic regression diagnostics. The Annals of Statistics, 9, 705-724.

Wei B.C. (1998). Exponential Family Nonlinear Models. Springer, Singapore.

Examples

###### Example 1: The effects of fertilizers on coastal Bermuda grass
data(Grass)
fit1 <- gnm(Yield ~ b0 + b1/(Nitrogen + a1) + b2/(Phosphorus + a2) + b3/(Potassium + a3),
            family=gaussian(inverse), start=c(b0=0.1,b1=13,b2=1,b3=1,a1=45,a2=15,a3=30), data=Grass)

fit1a <- update(fit1, subset=-c(1), start=coef(fit1), maxit=1)
coef(fit1) - coef(fit1a)

dfbetas <- dfbeta(fit1)
round(dfbetas[1,],5)

###### Example 2: Assay of an Insecticide with a Synergist
data(Melanopus)
fit2 <- gnm(Killed/Exposed ~ b0 + b1*log(Insecticide-a1) + b2*Synergist/(a2 + Synergist),
            family=binomial(logit), weights=Exposed, start=c(b0=-3,b1=1.2,a1=1.7,b2=1.7,a2=2),
		   data=Melanopus)

fit2a <- update(fit2, subset=-c(2), start=coef(fit2), maxit=1)
coef(fit2) - coef(fit2a)

dfbetas <- dfbeta(fit2)
round(dfbetas[2,],5)

###### Example 3: Developmental rate of Drosophila melanogaster
data(Drosophila)
fit3 <- gnm(Duration ~ b0 + b1*Temp + b2/(Temp-a), family=Gamma(log),
            start=c(b0=3,b1=-0.25,b2=-210,a=55), weights=Size, data=Drosophila)

fit3a <- update(fit3, subset=-c(3), start=coef(fit3), maxit=1)
coef(fit3) - coef(fit3a)

dfbetas <- dfbeta(fit3)
round(dfbetas[3,],5)

###### Example 4: Radioimmunological Assay of Cortisol
data(Cortisol)
fit4 <- gnm(Y ~ b0 + (b1-b0)/(1 + exp(b2+ b3*lDose))^b4, family=Gamma(identity),
            start=c(b0=130,b1=2800,b2=3,b3=3,b4=0.5), data=Cortisol)

fit4a <- update(fit4, subset=-c(4), start=coef(fit4), maxit=1)
coef(fit4) - coef(fit4a)

dfbetas <- dfbeta(fit4)
round(dfbetas[4,],5)

Description

Produces an approximation, better known as the one-step approximation, of the effect on the parameter estimates of deleting each individual in turn. This function also can produce an index plot of the Dfbeta statistic for some parameter chosen via the argument coefs.

Usage

## S3 method for class 'overglm'
dfbeta(model, coefs, identify, ...)

Arguments

Details

The one-step approximation of the estimates of the parameters when the i-th individual is excluded from the dataset consists of the vector obtained as result of the first iteration of the Newthon-Raphson algorithm when it is performed using: (1) a dataset in which the i-th individual is excluded; and (2) a starting value which is the estimate of the same model but based on the dataset inluding all individuals.

Value

A matrix with so many rows as individuals in the sample and so many columns as parameters in the linear predictor. The $i$-th row of that matrix corresponds to the difference between the estimates of the parameters in the linear predictor using all individuals and the one-step approximation of those estimates when the i-th individual is excluded from the dataset.

References

Pregibon D. (1981). Logistic regression diagnostics. The Annals of Statistics, 9, 705-724.

Examples

###### Example 1: Self diagnozed ear infections in swimmers
data(swimmers)
fit1 <- overglm(infections ~ frequency + location, family="nb1(log)", data=swimmers)
dfbeta(fit1, coefs="frequency", col="red", lty=1, lwd=1, col.lab="blue",
       col.axis="blue", col.main="black", family="mono", cex=0.8, main="frequency")

###### Example 2: Article production by graduate students in biochemistry PhD programs
bioChemists <- pscl::bioChemists
fit2 <- overglm(art ~ fem + kid5 + ment, family="nb1(log)", data = bioChemists)
dfbeta(fit2, coefs="fem", col="red", lty=1, lwd=1, col.lab="blue",
       col.axis="blue", col.main="black", family="mono", cex=0.8, main="fem")

###### Example 3: Agents to stimulate cellular differentiation
data(cellular)
fit3 <- overglm(cbind(cells,200-cells) ~ tnf + ifn, family="bb(logit)", data=cellular)
dfbeta(fit3, coefs="tnf", col="red", lty=1, lwd=1, col.lab="blue",
       col.axis="blue", col.main="black", family="mono", cex=0.8, main="tnf")

Description

Produces an approximation, better known as the one-step approximation, of the effect on the parameter estimates of deleting each individual in turn. This function also can produce an index plot of the Dfbeta statistic for some parameter chosen via the argument coefs.

Usage

## S3 method for class 'zeroinflation'
dfbeta(model, submodel = c("counts", "zeros"), coefs, identify, ...)

Arguments

Details

The one-step approximation of the estimates of the parameters when the i-th individual is excluded from the dataset consists of the vector obtained as result of the first iteration of the Newthon-Raphson algorithm when it is performed using: (1) a dataset in which the i-th individual is excluded; and (2) a starting value which is the estimate of the same model but based on the dataset inluding all individuals.

Value

A matrix with so many rows as individuals in the sample and so many columns as parameters in the linear predictor. The $i$-th row of that matrix corresponds to the difference between the estimates of the parameters in the linear predictor using all individuals and the one-step approximation of those estimates when the i-th individual is excluded from the dataset.

References

Pregibon D. (1981). Logistic regression diagnostics. The Annals of Statistics, 9, 705-724.

Examples

####### Example 1: Self diagnozed ear infections in swimmers
data(swimmers)
fit <- zeroinf(infections ~ frequency + location, family="nb1(log)", data=swimmers)

dfbeta(fit, submodel="counts", coefs="frequency", col="red", lty=1, lwd=1,
       col.lab="blue", col.axis="blue", col.main="black", family="mono", cex=0.8)

dfbeta(fit, submodel="zeros", coefs="location", col="red", lty=1, lwd=1,
       col.lab="blue", col.axis="blue", col.main="black", family="mono", cex=0.8)

Description

These data are counts of virus particles at 5 different dilutions. There are 4 replicate counts at each dilution except the last for which there are 5 counts. The aim is to estimate the expected number of virus particles per unit volume.

Usage

data(dilution)

Format

A data frame with 21 rows and 2 variables:

Count

a numeric vector indicating the count of virus particles.

Dilution

a numeric vector indicating the dilution volume.

Source

https://sada2013.sciencesconf.org/16138/glmSession4_Cotonou.pdf

Examples

data(dilution)
xlab <- "Dilution volume"
ylab <- "Count of virus particles"
dev.new()
with(dilution,plot(Dilution,Count,pch=20,xlab=xlab,ylab=ylab))

Description

Drosophila melanogaster developmental stages were monitored as part of an experiment to determine the effect of temperature on their duration. The eggs were laid at approximately 25 degrees Celsius and remained at that temperature for 20-30 minutes. The eggs were then brought to the experimental temperature, which remained constant throughout the experiment.

Usage

data(Drosophila)

Format

A data frame with 23 rows and 3 variables:

Temp

a numeric vector indicating the temperature, in degrees Celsius.

Duration

a numeric vector indicating the average duration of the embryonic period, in hours, measured from the time at which the eggs were laid.

Size

a numeric vector indicating how many eggs each batch contained.

References

Powsner L. (1935) The effects of temperature on the durations of the developmental stages of Drosophila melanogaster. Physiological Zoology, 8, 474-520.

McCullagh P., Nelder J.A. (1989). Generalized Linear Models. 2nd Edition. Chapman and Hall, London.

Wei B.C. (1998). Exponential Family Nonlinear Models. Springer, Singapore.

Examples

data(Drosophila)
dev.new()
with(Drosophila, plot(Temp,Duration,xlab="Temperature, in degrees Celsius",pch=16,col="blue",
                      ylab="Average duration of the embryonic period"))

Description

Generic function for building a normal QQ-plot with simulated envelope of residuals obtained from a fitted model.

Usage

envelope(object, ...)

Arguments

Value

A matrix with the simulated envelope and, optionally, a plot of it.

Description

Produces a normal QQ-plot with simulated envelope of residuals for generalized linear models.

Usage

## S3 method for class 'glm'
envelope(
  object,
  rep = 25,
  conf = 0.95,
  type = c("quantile", "deviance", "pearson"),
  standardized = FALSE,
  plot.it = TRUE,
  identify,
  ...
)

Arguments

Details

In order to construct the simulated envelope, rep independent realizations of the response variable for each individual are simulated, which is done by considering (1) the model assumption about the distribution of the response variable; (2) the estimation of the linear predictor parameters; and (3) the estimation of the dispersion parameter. Each time, the vector of observed responses is replaced with one of the simulated samples, re-fitting the interest model rep times. For each $i=1,2,...,n$, where $n$ is the number of individuals in the sample, the $i$-th order statistic of the type-type residuals is computed and then sorted for each replicate, giving a random sample of size rep of the $i$-th order statistic. In other words, the simulated envelope is comprised of the quantiles (1 - conf)/2 and (1 + conf)/2 of the random sample of size rep of the $i$-th order statistic of the type-type residuals for $i=1,2,...,n$.

Value

A matrix with the following four columns:

References

Atkinson A.C. (1985) Plots, Transformations and Regression. Oxford University Press, Oxford.

Davison A.C., Gigli A. (1989) Deviance Residuals and Normal Scores Plots. Biometrika 76, 211-221.

Dunn P.K., Smyth G.K. (1996) Randomized Quantile Residuals. Journal of Computational and Graphical Statistics 5, 236-244.

Pierce D.A., Schafer D.W. (1986) Residuals in Generalized Linear Models. Journal of the American Statistical Association 81, 977-986.

Wei B.C. (1998). Exponential Family Nonlinear Models. Springer, Singapore.

See Also

envelope.lm, envelope.gnm, envelope.overglm

Examples

###### Example 1:
burn1000 <- aplore3::burn1000
burn1000 <- within(burn1000, death <- factor(death, levels=c("Dead","Alive")))
fit1 <- glm(death ~ age*inh_inj + tbsa*inh_inj, family=binomial("logit"), data=burn1000)
envelope(fit1, rep=50, conf=0.95, type="pearson", col="red", pch=20, col.lab="blue",
         col.axis="blue", col.main="black", family="mono", cex=0.8)

###### Example 2: Fuel consumption of automobiles
Auto <- ISLR::Auto
fit2 <- glm(mpg ~ horsepower*weight, family=inverse.gaussian("log"), data=Auto)
envelope(fit2, rep=50, conf=0.95, type="pearson", col="red", pch=20, col.lab="blue",
         col.axis="blue", col.main="black", family="mono", cex=0.8)

###### Example 3: Skin cancer in women
data(skincancer)
fit3 <- glm(cases ~ city + ageC, offset=log(population), family=poisson, data=skincancer)
envelope(fit3, rep=100, conf=0.95, type="quantile", col="red", pch=20,col.lab="blue",
         col.axis="blue",col.main="black",family="mono",cex=0.8)

###### Example 4: Self diagnozed ear infections in swimmers
data(swimmers)
fit4 <- glm(infections ~ frequency + location, family=poisson(log), data=swimmers)
envelope(fit4, rep=100, conf=0.95, type="quantile", col="red", pch=20, col.lab="blue",
         col.axis="blue", col.main="black", family="mono", cex=0.8)

###### Example 5: Agents to stimulate cellular differentiation
data(cellular)
fit5 <- glm(cbind(cells,200-cells) ~ tnf + ifn, family=binomial(logit), data=cellular)
envelope(fit5, rep=100, conf=0.95, type="quantile", col="red", pch=20, col.lab="blue",
         col.axis="blue", col.main="black", family="mono", cex=0.8)

Description

Produces a normal QQ-plot with simulated envelope of residuals for generalized nonlinear models.

Usage

## S3 method for class 'gnm'
envelope(
  object,
  rep = 25,
  conf = 0.95,
  type = c("quantile", "deviance", "pearson"),
  standardized = FALSE,
  plot.it = TRUE,
  identify,
  ...
)

Arguments

Details

In order to construct the simulated envelope, rep independent realizations of the response variable for each individual are simulated, which is done by considering (1) the model assumption about the distribution of the response variable; (2) the estimation of the "linear" predictor parameters; and (3) the estimation of the dispersion parameter. Each time, the vector of observed responses is replaced with one of the simulated samples, re-fitting the interest model rep times. For each $i=1,2,...,n$, where $n$ is the number of individuals in the sample, the $i$-th order statistic of the type-type residuals is computed and then sorted for each replicate, giving a random sample of size rep of the $i$-th order statistic. In other words, the simulated envelope is comprised of the quantiles (1 - conf)/2 and (1 + conf)/2 of the random sample of size rep of the $i$-th order statistic of the type-type residuals for $i=1,2,...,n$.