library(contrast)
## Warning: package 'contrast' was built under R version 3.4.4
## Loading required package: rms
## Warning: package 'rms' was built under R version 3.4.4
## Loading required package: Hmisc
## Warning: package 'Hmisc' was built under R version 3.4.4
## Loading required package: lattice
## Loading required package: survival
## Loading required package: Formula
## Loading required package: ggplot2
##
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:base':
##
## format.pval, units
## Loading required package: SparseM
##
## Attaching package: 'SparseM'
## The following object is masked from 'package:base':
##
## backsolve
library(multcomp)
## Warning: package 'multcomp' was built under R version 3.4.4
## Loading required package: mvtnorm
## Loading required package: TH.data
## Warning: package 'TH.data' was built under R version 3.4.2
## Loading required package: MASS
## Warning: package 'MASS' was built under R version 3.4.3
##
## Attaching package: 'TH.data'
## The following object is masked from 'package:MASS':
##
## geyser
library(VGAM)
## Warning: package 'VGAM' was built under R version 3.4.4
## Loading required package: stats4
## Loading required package: splines
##
## Attaching package: 'VGAM'
## The following objects are masked from 'package:rms':
##
## calibrate, lrtest
Read-in Data
contra <- read.csv("C:/Users/miche/Documents/BS 853 Generalized Linear Model/contraceptive.csv", header= T)
dim(contra)
## [1] 21 5
Data Structure
#Nominal:1=Sterilization,2=Other,3=None#
contra$cage=contra$age-4;
contra$agesq=contra$cage^2;
contra$elsn=log((contra$count/contra$total)/(contra$ref/contra$total));
Postscript Set Up
#Encapsulated PostScript#
setEPS(paper='special')
postscript('contraceptive.eps',width=7,height=5)
plot(contra$age[contra$method %in% c(1,2)], contra$elsn[contra$method %in% c(1,2)],xlab='AGE',ylab='Log Generalized Logit (None as reference)',type='n');
lines(contra$age[contra$method %in% c(1)], contra$elsn[contra$method %in% c(1)],cex=1.4,type='b');
lines(contra$age[contra$method %in% c(2)], contra$elsn[contra$method %in% c(2)],cex=1.4,type='b',pch=19);
legend('bottomright',legend=c('Sterilization', 'Other'), title='Method',pch=c(1,19),pt.cex=1.4)
dev.off()
## png
## 2
Factors
contra$method <- factor(contra$method)
contra$age<- factor(contra$age,levels=c(7,1:6))
Model 1 Using vglm()
cat('\n Contraceptive data \n Using MLogit - AGE categorical \n');
##
## Contraceptive data
## Using MLogit - AGE categorical
#vglm function fits wide variety of models.
#These models include the cumulative logit model with proportional odds.
M<-vglm(method ~ age, multinomial, weight=count, contra)
summary(M)
##
## Call:
## vglm(formula = method ~ age, family = multinomial, data = contra,
## weights = count)
##
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## log(mu[,1]/mu[,3]) -15.78 -9.107 -2.030 13.88 23.02
## log(mu[,2]/mu[,3]) -11.78 -5.576 -2.002 16.48 22.30
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 -0.6986 0.1283 -5.446 5.14e-08 ***
## (Intercept):2 -2.9069 0.3248 -8.951 < 2e-16 ***
## age1:1 -3.6495 0.5941 -6.142 8.12e-10 ***
## age1:2 1.5710 0.3552 4.423 9.74e-06 ***
## age2:1 -0.9108 0.1774 -5.136 2.81e-07 ***
## age2:2 1.8354 0.3395 5.406 6.44e-08 ***
## age3:1 0.3668 0.1562 2.348 0.018879 *
## age3:2 2.0750 0.3412 6.081 1.19e-09 ***
## age4:1 0.9764 0.1585 6.162 7.20e-10 ***
## age4:2 1.9244 0.3515 5.475 4.38e-08 ***
## age5:1 0.7454 0.1639 4.549 5.39e-06 ***
## age5:2 1.5825 0.3616 4.376 1.21e-05 ***
## age6:1 0.6094 0.1709 3.565 0.000364 ***
## age6:2 0.9851 0.3914 2.517 0.011851 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Number of linear predictors: 2
##
## Names of linear predictors: log(mu[,1]/mu[,3]), log(mu[,2]/mu[,3])
##
## Residual deviance: 5745.798 on 28 degrees of freedom
##
## Log-likelihood: -2872.899 on 28 degrees of freedom
##
## Number of iterations: 6
##
## Warning: Hauck-Donner effect detected in the following estimate(s):
## '(Intercept):2', 'age1:1'
##
## Reference group is level 3 of the response
LogLinear Model
cat('\n Contraceptive data \n Using Loglinear models - AGE categorical \n');
##
## Contraceptive data
## Using Loglinear models - AGE categorical
M2<- glm(count ~ method*age, family=poisson(), data=contra);
Model 2
cat('\n Contraceptive data \n Using MLogit - - continuous age linear effect \n');
##
## Contraceptive data
## Using MLogit - - continuous age linear effect
M<-vglm(method~cage,data=contra,multinomial,weights=count)
summary(M)
##
## Call:
## vglm(formula = method ~ cage, family = multinomial, data = contra,
## weights = count)
##
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## log(mu[,1]/mu[,3]) -15.30 -10.874 -3.001 8.754 24.48
## log(mu[,2]/mu[,3]) -12.19 -6.344 -1.892 13.059 24.62
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 -0.49983 0.04098 -12.196 < 2e-16 ***
## (Intercept):2 -1.35799 0.05847 -23.225 < 2e-16 ***
## cage:1 0.26714 0.02340 11.415 < 2e-16 ***
## cage:2 -0.18556 0.03230 -5.745 9.18e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Number of linear predictors: 2
##
## Names of linear predictors: log(mu[,1]/mu[,3]), log(mu[,2]/mu[,3])
##
## Residual deviance: 6039.776 on 38 degrees of freedom
##
## Log-likelihood: -3019.888 on 38 degrees of freedom
##
## Number of iterations: 5
##
## No Hauck-Donner effect found in any of the estimates
##
## Reference group is level 3 of the response
Model 3
cat('\n Contraceptive data \n Using MLogit - continuous age linear and quadratic effects\n');
##
## Contraceptive data
## Using MLogit - continuous age linear and quadratic effects
M<-vglm(method~cage+agesq,data=contra,multinomial,weights=count)
summary(M)
##
## Call:
## vglm(formula = method ~ cage + agesq, family = multinomial, data = contra,
## weights = count)
##
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## log(mu[,1]/mu[,3]) -15.46 -8.666 -2.999 13.30 22.80
## log(mu[,2]/mu[,3]) -11.89 -5.626 -1.780 15.93 22.38
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 0.16749 0.05817 2.879 0.00399 **
## (Intercept):2 -0.99314 0.07634 -13.010 < 2e-16 ***
## cage:1 0.38545 0.02950 13.066 < 2e-16 ***
## cage:2 -0.22600 0.03705 -6.100 1.06e-09 ***
## agesq:1 -0.24332 0.01647 -14.773 < 2e-16 ***
## agesq:2 -0.11895 0.01899 -6.264 3.76e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Number of linear predictors: 2
##
## Names of linear predictors: log(mu[,1]/mu[,3]), log(mu[,2]/mu[,3])
##
## Residual deviance: 5766.273 on 36 degrees of freedom
##
## Log-likelihood: -2883.136 on 36 degrees of freedom
##
## Number of iterations: 5
##
## No Hauck-Donner effect found in any of the estimates
##
## Reference group is level 3 of the response
Using LogLinear Model
cat('\n Contraceptive data \n Using LogLinear models - continuous age linear and quadratic effects\n');
##
## Contraceptive data
## Using LogLinear models - continuous age linear and quadratic effects
M2<- glm(count ~ method+age+method:cage + method:agesq, family=poisson(), data=contra)
summary(M2)
##
## Call:
## glm(formula = count ~ method + age + method:cage + method:agesq,
## family = poisson(), data = contra)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.47791 -0.59984 0.01931 0.48106 1.69573
##
## Coefficients: (2 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.42099 0.09486 57.149 < 2e-16 ***
## method2 -1.16063 0.07958 -14.585 < 2e-16 ***
## method3 -0.16749 0.05817 -2.879 0.00399 **
## age1 0.18083 0.09185 1.969 0.04898 *
## age2 0.72134 0.08243 8.751 < 2e-16 ***
## age3 0.50610 0.08558 5.914 3.34e-09 ***
## age4 0.11378 0.08975 1.268 0.20488
## age5 -0.14334 0.08823 -1.625 0.10426
## age6 -0.17818 0.08343 -2.136 0.03271 *
## method1:cage 0.38545 0.02950 13.066 < 2e-16 ***
## method2:cage -0.22600 0.03705 -6.099 1.07e-09 ***
## method3:cage NA NA NA NA
## method1:agesq -0.24332 0.01647 -14.773 < 2e-16 ***
## method2:agesq -0.11895 0.01899 -6.264 3.76e-10 ***
## method3:agesq NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 1521.567 on 20 degrees of freedom
## Residual deviance: 20.475 on 8 degrees of freedom
## AIC: 182.06
##
## Number of Fisher Scoring iterations: 4