gnm
Generalized Nonlinear Models

Package index
Search the gnm package
Vignettes
Functions
125
Source code
129
Man pages
48
Home
/
CRAN
/
gnm
/
inst/doc/gnmOverview.R

inst/doc/gnmOverview.R
In gnm: Generalized Nonlinear Models 

### R code from vignette source 'gnmOverview.Rnw'

###################################################
### code chunk number 1: gnmOverview.Rnw:65-66
###################################################
getOption("SweaveHooks")[["eval"]]()
options(SweaveHooks = list(eval = function() options(show.signif.stars = FALSE)))


###################################################
### code chunk number 2: Load_gnm
###################################################
getOption("SweaveHooks")[["eval"]]()
library(gnm)


###################################################
### code chunk number 3: migrationData
###################################################
getOption("SweaveHooks")[["eval"]]()
count <- c(11607,   100,   366,   124,
              87, 13677,   515,   302,
             172,   225, 17819,   270,
              63,   176,   286, 10192 )
region <- c("NE", "MW", "S", "W")
row <-  gl(4, 4, labels = region)
col <-  gl(4, 1, length = 16, labels = region)


###################################################
### code chunk number 4: squareTableModels
###################################################
getOption("SweaveHooks")[["eval"]]()
independence <- glm(count ~ row + col, family = poisson)
quasi.indep <- glm(count ~ row + col + Diag(row, col), family = poisson)
symmetry <- glm(count ~ Symm(row, col), family = poisson)
quasi.symm <- glm(count ~ row + col + Symm(row, col), family = poisson)
comparison1 <- anova(independence, quasi.indep, quasi.symm)
print(comparison1, digits = 7)
comparison2 <- anova(symmetry, quasi.symm)
print(comparison2)


###################################################
### code chunk number 5: EriksonData
###################################################
getOption("SweaveHooks")[["eval"]]()
### Collapse to 7 by 7 table as in Erikson et al. (1982)
erikson <- as.data.frame(erikson)
lvl <- levels(erikson$origin)
levels(erikson$origin) <- levels(erikson$destination) <-
    c(rep(paste(lvl[1:2], collapse = " + "), 2), lvl[3],
      rep(paste(lvl[4:5], collapse = " + "), 2), lvl[6:9])
erikson <- xtabs(Freq ~ origin + destination + country, data = erikson)


###################################################
### code chunk number 6: wedderburn
###################################################
getOption("SweaveHooks")[["eval"]]()
##  data from Wedderburn (1974), see ?barley
logitModel <- glm(y ~ site + variety, family = wedderburn, data = barley)
fit <- fitted(logitModel)
print(sum((barley$y - fit)^2 / (fit * (1-fit))^2))


###################################################
### code chunk number 7: termPredictors
###################################################
getOption("SweaveHooks")[["eval"]]()
print(temp <- termPredictors(quasi.symm))
rowSums(temp) - quasi.symm$linear.predictors


###################################################
### code chunk number 8: RC_homogeneous_model_1
###################################################
getOption("SweaveHooks")[["eval"]]()
set.seed(1)
RChomog1 <- gnm(Freq ~ origin + destination + Diag(origin, destination) +
               MultHomog(origin, destination), family = poisson,
               data = occupationalStatus, verbose = FALSE)


###################################################
### code chunk number 9: RC_homogeneous_model_2
###################################################
getOption("SweaveHooks")[["eval"]]()
set.seed(2)
RChomog2 <- update(RChomog1)


###################################################
### code chunk number 10: Compare_coefficients
###################################################
getOption("SweaveHooks")[["eval"]]()
compareCoef <- cbind(coef(RChomog1), coef(RChomog2))
colnames(compareCoef) <- c("RChomog1", "RChomog2")
round(compareCoef, 4)


###################################################
### code chunk number 11: Summarize_model
###################################################
getOption("SweaveHooks")[["eval"]]()
summary(RChomog2)


###################################################
### code chunk number 12: RC_homogeneous_constrained_model1
###################################################
getOption("SweaveHooks")[["eval"]]()
set.seed(1)
RChomogConstrained1 <- update(RChomog1, constrain = length(coef(RChomog1)))


###################################################
### code chunk number 13: RC_homogeneous_constrained_model2
###################################################
getOption("SweaveHooks")[["eval"]]()
set.seed(2)
RChomogConstrained2 <- update(RChomogConstrained1)
identical(coef(RChomogConstrained1), coef(RChomogConstrained2))


###################################################
### code chunk number 14: Eliminate_Eg
###################################################
getOption("SweaveHooks")[["eval"]]()
set.seed(1)
n <- 1000
x <- rep(rnorm(n), rep(3, n))
counts <- as.vector(rmultinom(n, 10, c(0.7, 0.1, 0.2)))
rowID <- gl(n, 3, 3 * n)
resp <- gl(3, 1, 3 * n)


###################################################
### code chunk number 15: Double_UNIDIFF_model
###################################################
getOption("SweaveHooks")[["eval"]]()
doubleUnidiff <- gnm(Freq ~ election:vote + election:class:religion
                     + Mult(Exp(election), religion:vote) +
                     Mult(Exp(election), class:vote), family = poisson,
                     data = cautres)


###################################################
### code chunk number 16: Contrast_matrix
###################################################
getOption("SweaveHooks")[["eval"]]()
coefs <- names(coef(doubleUnidiff))
contrCoefs <- coefs[grep(", religion:vote", coefs)]
nContr <- length(contrCoefs)
contrMatrix <- matrix(0, length(coefs), nContr,
                      dimnames = list(coefs, contrCoefs))
contr <- contr.sum(contrCoefs)
# switch round to contrast with first level
contr <- rbind(contr[nContr, ], contr[-nContr, ])
contrMatrix[contrCoefs, 2:nContr] <- contr
contrMatrix[contrCoefs, 2:nContr]


###################################################
### code chunk number 17: Check_estimability_1
###################################################
getOption("SweaveHooks")[["eval"]]()
checkEstimable(doubleUnidiff, contrMatrix)


###################################################
### code chunk number 18: Check_estimability_2
###################################################
getOption("SweaveHooks")[["eval"]]()
coefs <- names(coef(doubleUnidiff))
contrCoefs <- coefs[grep("[.]religion", coefs)]
nContr <- length(contrCoefs)
contrMatrix <- matrix(0, length(coefs), length(contrCoefs),
                      dimnames = list(coefs, contrCoefs))
contr <- contr.sum(contrCoefs)
contrMatrix[contrCoefs, 2:nContr] <- rbind(contr[nContr, ], contr[-nContr, ])
checkEstimable(doubleUnidiff, contrMatrix)


###################################################
### code chunk number 19: Get_contrasts_1
###################################################
getOption("SweaveHooks")[["eval"]]()
myContrasts <- getContrasts(doubleUnidiff,
                            pickCoef(doubleUnidiff, ", religion:vote"))
myContrasts


###################################################
### code chunk number 20: qvplot
###################################################
getOption("SweaveHooks")[["eval"]]()
plot(myContrasts,
  main = "Relative strength of religion-vote association, log scale",
xlab = "Election", levelNames = 1:4)


###################################################
### code chunk number 21: RCmodel
###################################################
getOption("SweaveHooks")[["eval"]]()
mentalHealth$MHS <- C(mentalHealth$MHS, treatment)
mentalHealth$SES <- C(mentalHealth$SES, treatment)
RC1model <- gnm(count ~ SES + MHS + Mult(SES, MHS),
                family = poisson, data = mentalHealth)


###################################################
### code chunk number 22: RCmodel_constrained
###################################################
getOption("SweaveHooks")[["eval"]]()
RC1model2 <- gnm(count ~ SES + MHS + Mult(1, SES, MHS),
                 constrain = "[.]SES[AF]", constrainTo = c(0, 1),
                 ofInterest = "[.]SES",
                 family = poisson, data = mentalHealth)
summary(RC1model2)


###################################################
### code chunk number 23: getContrasts_simple
###################################################
getOption("SweaveHooks")[["eval"]]()
getContrasts(RC1model, pickCoef(RC1model, "[.]SES"), ref = "first",
             scaleRef = "first", scaleWeights = c(rep(0, 5), 1))


###################################################
### code chunk number 24: two-way
###################################################
getOption("SweaveHooks")[["eval"]]()
xtabs(y ~ site + variety, barley)


###################################################
### code chunk number 25: residSVD
###################################################
getOption("SweaveHooks")[["eval"]]()
emptyModel <- gnm(y ~ -1, family = wedderburn, data = barley)
biplotStart <- residSVD(emptyModel, barley$site, barley$variety, d = 2)
biplotModel <- gnm(y ~ -1 + instances(Mult(site, variety), 2),
family = wedderburn, data = barley, start =  biplotStart)


###################################################
### code chunk number 26: residSVDplot
###################################################
getOption("SweaveHooks")[["eval"]]()
plot(coef(biplotModel), biplotStart,
     main = "Comparison of residSVD and MLE for a 2-dimensional
 biplot model", ylim = c(-2, 2), xlim = c(-4, 4))
abline(a = 0, b = 1, lty = 2)


###################################################
### code chunk number 27: Set_contrasts_attribute
###################################################
getOption("SweaveHooks")[["eval"]]()
set.seed(1)
mentalHealth$MHS <- C(mentalHealth$MHS, treatment)
mentalHealth$SES <- C(mentalHealth$SES, treatment)


###################################################
### code chunk number 28: RC1_model
###################################################
getOption("SweaveHooks")[["eval"]]()
RC1model <- gnm(count ~ SES + MHS + Mult(SES, MHS), family = poisson,
                data = mentalHealth)
RC1model


###################################################
### code chunk number 29: Normalize_scores
###################################################
getOption("SweaveHooks")[["eval"]]()
rowProbs <- with(mentalHealth, tapply(count, SES, sum) / sum(count))
colProbs <- with(mentalHealth, tapply(count, MHS, sum) / sum(count))
rowScores <- coef(RC1model)[10:15]
colScores <- coef(RC1model)[16:19]
rowScores <- rowScores - sum(rowScores * rowProbs)
colScores <- colScores - sum(colScores * colProbs)
beta1 <- sqrt(sum(rowScores^2 * rowProbs))
beta2 <- sqrt(sum(colScores^2 * colProbs))
assoc <- list(beta = beta1 * beta2,
              mu = rowScores / beta1,
              nu = colScores / beta2)
assoc


###################################################
### code chunk number 30: Elliptical_contrasts
###################################################
getOption("SweaveHooks")[["eval"]]()
mu <- getContrasts(RC1model, pickCoef(RC1model, "[.]SES"),
                   ref = rowProbs, scaleWeights = rowProbs)
nu <- getContrasts(RC1model, pickCoef(RC1model, "[.]MHS"),
                   ref = colProbs, scaleWeights = colProbs)
mu
nu


###################################################
### code chunk number 31: RC2_model
###################################################
getOption("SweaveHooks")[["eval"]]()
RC2model <- gnm(count ~ SES + MHS + instances(Mult(SES, MHS), 2),
                family = poisson, data = mentalHealth)
RC2model


###################################################
### code chunk number 32: Homogeneous_effects
###################################################
getOption("SweaveHooks")[["eval"]]()
RChomog <- gnm(Freq ~ origin + destination + Diag(origin, destination) +
               MultHomog(origin, destination), family = poisson,
               data = occupationalStatus)
RChomog


###################################################
### code chunk number 33: Heterogeneous_effects
###################################################
getOption("SweaveHooks")[["eval"]]()
RCheterog <- gnm(Freq ~ origin + destination + Diag(origin, destination) +
               Mult(origin, destination), family = poisson,
               data = occupationalStatus)
anova(RChomog, RCheterog)


###################################################
### code chunk number 34: Transform_to_counts
###################################################
getOption("SweaveHooks")[["eval"]]()
set.seed(1)
count <- with(voting, percentage/100 * total)
yvar <- cbind(count, voting$total - count)


###################################################
### code chunk number 35: Class_mobility
###################################################
getOption("SweaveHooks")[["eval"]]()
classMobility <- gnm(yvar ~ Dref(origin, destination),
                       family = binomial, data = voting)
classMobility


###################################################
### code chunk number 36: Class_mobility_weights
###################################################
getOption("SweaveHooks")[["eval"]]()
DrefWeights(classMobility)


###################################################
### code chunk number 37: Salariat_factors
###################################################
getOption("SweaveHooks")[["eval"]]()
upward <- with(voting, origin != 1 & destination == 1)
downward <- with(voting, origin == 1 & destination != 1)


###################################################
### code chunk number 38: Social_mobility
###################################################
getOption("SweaveHooks")[["eval"]]()
socialMobility <- gnm(yvar ~ Dref(origin, destination,
                                  delta = ~ 1 + downward + upward),
                      family = binomial, data = voting)
socialMobility


###################################################
### code chunk number 39: social_mobility_weights
###################################################
getOption("SweaveHooks")[["eval"]]()
DrefWeights(socialMobility)


###################################################
### code chunk number 40: Downward_mobility
###################################################
getOption("SweaveHooks")[["eval"]]()
downwardMobility <- gnm(yvar ~ Dref(origin, destination,
                                    delta = ~ 1 + downward),
                        family = binomial, data = voting)
downwardMobility
DrefWeights(downwardMobility)


###################################################
### code chunk number 41: UNIDIFF_model
###################################################
getOption("SweaveHooks")[["eval"]]()
set.seed(1)
unidiff <- gnm(Freq ~ educ*orig + educ*dest + Mult(Exp(educ), orig:dest),
               ofInterest = "[.]educ", family = poisson,
               data = yaish, subset = (dest != 7))
coef(unidiff)


###################################################
### code chunk number 42: Unidiff_contrasts
###################################################
getOption("SweaveHooks")[["eval"]]()
getContrasts(unidiff, ofInterest(unidiff))


###################################################
### code chunk number 43: double_UNIDIFF_model
###################################################
getOption("SweaveHooks")[["eval"]]()
set.seed(1)
doubleUnidiff <- gnm(Freq ~ election*vote + election*class*religion +
                     Mult(Exp(election), religion:vote) +
                     Mult(Exp(election), class:vote),

             family = poisson, data = cautres)
getContrasts(doubleUnidiff, rev(pickCoef(doubleUnidiff, ", class:vote")))
getContrasts(doubleUnidiff, rev(pickCoef(doubleUnidiff, ", religion:vote")))


###################################################
### code chunk number 44: Scale_yields
###################################################
getOption("SweaveHooks")[["eval"]]()
set.seed(1)
yield.scaled <- wheat$yield * sqrt(3/1000)
treatment <- interaction(wheat$tillage, wheat$summerCrop, wheat$manure,
                         wheat$N, sep = "")


###################################################
### code chunk number 45: AMMI_model
###################################################
getOption("SweaveHooks")[["eval"]]()
mainEffects <- gnm(yield.scaled ~ year + treatment, family = gaussian,
                   data = wheat)
svdStart <- residSVD(mainEffects, year, treatment, 3)
bilinear1 <- update(mainEffects, . ~ . + Mult(year, treatment),
                    start = c(coef(mainEffects), svdStart[,1]))


###################################################
### code chunk number 46: AOD
###################################################
getOption("SweaveHooks")[["eval"]]()
anova(mainEffects, bilinear1, test = "F")


###################################################
### code chunk number 47: AMMI_model2
###################################################
getOption("SweaveHooks")[["eval"]]()
set.seed(1)
barleyModel <- gnm(height ~ year + genotype + Mult(year, genotype),
                   data = barleyHeights)


###################################################
### code chunk number 48: Spherical_contrasts
###################################################
getOption("SweaveHooks")[["eval"]]()
gamma <- getContrasts(barleyModel, pickCoef(barleyModel, "[.]y"),
                      ref = "mean", scaleWeights = "unit")
delta <- getContrasts(barleyModel, pickCoef(barleyModel, "[.]g"),
                      ref = "mean", scaleWeights = "unit")
gamma
delta


###################################################
### code chunk number 49: CI
###################################################
getOption("SweaveHooks")[["eval"]]()
gamma[[2]][,1] + (gamma[[2]][,2]) %o% c(-1.96, 1.96)
delta[[2]][,1] + (delta[[2]][,2]) %o% c(-1.96, 1.96)


###################################################
### code chunk number 50: SVD
###################################################
getOption("SweaveHooks")[["eval"]]()
svd(termPredictors(barleyModel)[, "Mult(year, genotype)"])$d


###################################################
### code chunk number 51: Biplot_model
###################################################
getOption("SweaveHooks")[["eval"]]()
set.seed(83)
biplotModel <- gnm(y ~ -1 + instances(Mult(site, variety), 2),
                   family = wedderburn, data = barley)


###################################################
### code chunk number 52: Row_and_column_scores
###################################################
getOption("SweaveHooks")[["eval"]]()
barleyMatrix <- xtabs(biplotModel$predictors ~ site + variety,
                      data = barley)
barleySVD <- svd(barleyMatrix)
A <- sweep(barleySVD$u, 2, sqrt(barleySVD$d), "*")[, 1:2]
B <- sweep(barleySVD$v, 2, sqrt(barleySVD$d), "*")[, 1:2]
rownames(A) <- levels(barley$site)
rownames(B) <- levels(barley$variety)
colnames(A) <- colnames(B) <- paste("Component", 1:2)
A
B


###################################################
### code chunk number 53: Biplot1
###################################################
getOption("SweaveHooks")[["eval"]]()
barleyCol <- c("red", "blue")
plot(rbind(A, B), pch = c(levels(barley$site), levels(barley$variety)),
     col = rep(barleyCol, c(nlevels(barley$site), nlevels(barley$variety))),
     xlim = c(-4, 4), ylim = c(-4, 4), main = "Biplot for barley data",
     xlab = "Component 1", ylab = "Component 2")
text(c(-3.5, -3.5), c(3.9, 3.6), c("sites: A-I","varieties: 1-9, X"),
     col = barleyCol, adj = 0)


###################################################
### code chunk number 54: Biplot2
###################################################
getOption("SweaveHooks")[["eval"]]()
plot(rbind(A, B), pch = c(levels(barley$site), levels(barley$variety)),
     col = rep(barleyCol, c(nlevels(barley$site), nlevels(barley$variety))),
     xlim = c(-4, 4), ylim = c(-4, 4), main = "Biplot for barley data",
     xlab = "Component 1", ylab = "Component 2")
text(c(-3.5, -3.5), c(3.9, 3.6), c("sites: A-I","varieties: 1-9, X"),
     col = barleyCol, adj = 0)
abline(a = 0, b = tan(pi/3))
abline(a = 0, b = -tan(pi/6))
abline(a = 2.6, b = tan(pi/3), lty = 2)
abline(a = 4.5, b = tan(pi/3), lty = 2)
abline(a = 1.3, b = -tan(pi/6), lty = 2)
text(2.8, 3.9, "v-axis", font = 3)
text(3.8, -2.7, "h-axis", font = 3)


###################################################
### code chunk number 55: Double_additive
###################################################
getOption("SweaveHooks")[["eval"]]()
variety.binary <- factor(match(barley$variety, c(2,3,6), nomatch = 0) > 0,
                        labels = c("rest", "2,3,6"))
doubleAdditive <- gnm(y ~ variety + Mult(site, variety.binary),
                      family = wedderburn, data = barley)


###################################################
### code chunk number 56: Compare_chi-squared
###################################################
getOption("SweaveHooks")[["eval"]]()
biplotModChiSq <- sum(residuals(biplotModel, type = "pearson")^2)
doubleAddChiSq <- sum(residuals(doubleAdditive, type = "pearson")^2)
c(doubleAddChiSq - biplotModChiSq,
  doubleAdditive$df.residual - biplotModel$df.residual)


###################################################
### code chunk number 57: Re-express_data
###################################################
getOption("SweaveHooks")[["eval"]]()
set.seed(1)
subset(backPain, x1 == 1 & x2 == 1 & x3 == 1)
backPainLong <- expandCategorical(backPain, "pain")
head(backPainLong)


###################################################
### code chunk number 58: Stereotype_model
###################################################
getOption("SweaveHooks")[["eval"]]()
oneDimensional <- gnm(count ~ pain + Mult(pain, x1 + x2 + x3),
                      eliminate = id, family = "poisson", data = backPainLong)
oneDimensional


###################################################
### code chunk number 59: Qualitative_model
###################################################
getOption("SweaveHooks")[["eval"]]()
threeDimensional  <- gnm(count ~ pain + pain:(x1 + x2 + x3), eliminate = id,
                         family = "poisson", data = backPainLong)
threeDimensional


###################################################
### code chunk number 60: Calculate_log-likelihood
###################################################
getOption("SweaveHooks")[["eval"]]()
logLikMultinom <- function(model, size){
    object <- get(model)
    l <- sum(object$y * log(object$fitted/size))
    c(nParameters = object$rank - nlevels(object$eliminate), logLikelihood = l)
}
size <- tapply(backPainLong$count, backPainLong$id, sum)[backPainLong$id]
t(sapply(c("oneDimensional", "threeDimensional"), logLikMultinom, size))


###################################################
### code chunk number 61: Constrain_slopes
###################################################
getOption("SweaveHooks")[["eval"]]()
## before constraint
summary(oneDimensional)
oneDimensional <- gnm(count ~ pain + Mult(pain, offset(x1) + x2 + x3),
                      eliminate = id, family = "poisson", data = backPainLong)
## after constraint
summary(oneDimensional)


###################################################
### code chunk number 62: Get_slopes
###################################################
getOption("SweaveHooks")[["eval"]]()
getContrasts(oneDimensional, pickCoef(oneDimensional, "[.]pain"))


###################################################
### code chunk number 63: singleExp
###################################################
getOption("SweaveHooks")[["eval"]]()
x <- 1:100
y <- exp(- x / 10)
set.seed(1)
saved.fits <- list()
for (i in 1:100) saved.fits[[i]] <- gnm(y ~ Exp(1 + x), verbose = FALSE)
table(zapsmall(sapply(saved.fits, deviance)))


###################################################
### code chunk number 64: singleExp2
###################################################
getOption("SweaveHooks")[["eval"]]()
saved.fits[[2]]


###################################################
### code chunk number 65: doubleExp
###################################################
getOption("SweaveHooks")[["eval"]]()
x <- 1:100
y <- exp(- x / 10) + 2 * exp(- x / 50)
set.seed(1)
saved.fits <- list()
for (i in 1:100) {
    saved.fits[[i]] <- suppressWarnings(gnm(y ~ Exp(1 + x, inst = 1) +
                                            Exp(1 + x, inst = 2),
                                            verbose = FALSE))
}
table(round(unlist(sapply(saved.fits, deviance)), 4))


###################################################
### code chunk number 66: doubleExp2
###################################################
getOption("SweaveHooks")[["eval"]]()
singleExp <- gnm(y ~ Exp(1 + x), start = c(NA, NA, -0.1), verbose = FALSE)
singleExp
meanOnly <- gnm(y ~ 1, verbose = FALSE)
meanOnly
plot(x, y, main = "Two sub-optimal fits to a sum-of-exponentials curve")
lines(x, fitted(singleExp))
lines(x, fitted(meanOnly), lty = "dashed")


###################################################
### code chunk number 67: doubleExp3
###################################################
getOption("SweaveHooks")[["eval"]]()
gnm(y ~ instances(Exp(1 + x), 2), start = c(NA, NA, -0.1, NA, -0.1),
    verbose = FALSE)

Try the gnm package in your browser

Any scripts or data that you put into this service are public.

gnm documentation built on Sept. 16, 2023, 5:06 p.m.

R Package Documentation

rdrr.io home R language documentation Run R code online

Browse R Packages

CRAN packages Bioconductor packages R-Forge packages GitHub packages

We want your feedback!

Note that we can't provide technical support on individual packages. You should contact the package authors for that.
GitHub issue tracker
ian@mutexlabs.com
Personal blog

 
Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.

Close