Nothing
### 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)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.