Class 2 Generalized Linear Model Using R: Short Practice

CareMod<-  glm(count ~ lessca*death, family=poisson(), Caredat)  #contains all components of the model
names(CareMod)

##  [1] "coefficients"      "residuals"         "fitted.values"    
##  [4] "effects"           "R"                 "rank"             
##  [7] "qr"                "family"            "linear.predictors"
## [10] "deviance"          "aic"               "null.deviance"    
## [13] "iter"              "weights"           "prior.weights"    
## [16] "df.residual"       "df.null"           "y"                
## [19] "converged"         "boundary"          "model"            
## [22] "call"              "formula"           "terms"            
## [25] "data"              "offset"            "control"          
## [28] "method"            "contrasts"         "xlevels"

CareMod$family 

## 
## Family: poisson 
## Link function: log

CareMod$contrasts

## $lessca
## [1] "contr.treatment"
## 
## $death
## [1] "contr.treatment"

S1<- summary(CareMod)

S1<- summary(CareMod)
S1$coefficients

##                  Estimate Std. Error    z value     Pr(>|z|)
## (Intercept)     5.7557422 0.05625440 102.316312 0.000000e+00
## lesscay         0.1658362 0.07645601   2.169041 3.007959e-02
## deathy         -3.9639827 0.41210584  -9.618846 6.657360e-22
## lesscay:deathy  1.0381366 0.47171197   2.200785 2.775125e-02

Change the contrats so that the output crrespond to the SAS output

options(contrasts = c("contr.SAS", "contr.poly"))
options(contrasts = c("contr.SAS", "contr.poly")) 
  CareMod2<-  glm(count ~ lessca*death, family=poisson(),data=Caredat) 
  CareMod2$family 

## 
## Family: poisson 
## Link function: log

  CareMod2$contrasts  

## $lessca
## [1] "contr.SAS"
## 
## $death
## [1] "contr.SAS"

  S2<- summary(CareMod2) 
  S2$coefficients  

##                 Estimate Std. Error   z value     Pr(>|z|)
## (Intercept)     2.995732  0.2236068 13.397322 6.268444e-41
## lesscan        -1.203973  0.4654747 -2.586549 9.694252e-03
## deathn          2.925846  0.2295233 12.747488 3.219516e-37
## lesscan:deathn  1.038137  0.4717120  2.200785 2.775125e-02

l<-sapply(Caredat, function(x) length(levels(x)))
###lessca  death  count 
###     2      2      0
T <- array(0, l[-3], lapply(Caredat[, -3], levels))
T[data.matrix(Caredat[,-3])] <- Caredat$count  

LNf <- loglm(~lessca*death,T);LNf0<-loglm(~1:2,T)



LL2 <-  glm(count ~ lessca + death, family=poisson(),Caredat)  
summary(LL2)  # Get the coefficients

## 
## Call:
## glm(formula = count ~ lessca + death, family = poisson(), data = Caredat)
## 
## Deviance Residuals: 
##       1        2        3        4  
##  1.4234  -1.8425  -0.2941   0.3231  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.65962    0.19901  13.364  < 2e-16 ***
## lesscan     -0.19926    0.07517  -2.651  0.00803 ** 
## deathn       3.27714    0.19978  16.404  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 780.494  on 3  degrees of freedom
## Residual deviance:   5.612  on 1  degrees of freedom
## AIC: 35.465
## 
## Number of Fisher Scoring iterations: 4

#              Estimate Std. Error   z value     Pr(>|z|)
#(Intercept)  2.6596236 0.19901495 13.363939 9.822600e-41
#lesscan     -0.1992581 0.07516726 -2.650863 8.028651e-03
#deathn       3.2771447 0.19978089 16.403695 1.799518e-60

LNi <- loglm(~lessca+death,T);LNi0<-loglm(~1+2,T)
LNi

## Call:
## loglm(formula = ~lessca + death, data = T)
## 
## Statistics:
##                       X^2 df   P(> X^2)
## Likelihood Ratio 5.611979  1 0.01783811
## Pearson          5.255464  1 0.02187797