Examples for 'stats::family'

Family Objects for Models

Aliases: family binomial gaussian Gamma inverse.gaussian poisson quasi quasibinomial quasipoisson

Keywords: models

### ** Examples

require(utils) # for str

nf <- gaussian()  # Normal family
nf

Family: gaussian 
Link function: identity 

str(nf)

List of 11
 $ family    : chr "gaussian"
 $ link      : chr "identity"
 $ linkfun   :function (mu)  
 $ linkinv   :function (eta)  
 $ variance  :function (mu)  
 $ dev.resids:function (y, mu, wt)  
 $ aic       :function (y, n, mu, wt, dev)  
 $ mu.eta    :function (eta)  
 $ initialize:  expression({  n <- rep.int(1, nobs)  if (is.null(etastart) && is.null(start) && is.null(mustart) && ((family$link| __truncated__
 $ validmu   :function (mu)  
 $ valideta  :function (eta)  
 - attr(*, "class")= chr "family"

gf <- Gamma()
gf

Family: Gamma 
Link function: inverse 

str(gf)

List of 12
 $ family    : chr "Gamma"
 $ link      : chr "inverse"
 $ linkfun   :function (mu)  
 $ linkinv   :function (eta)  
 $ variance  :function (mu)  
 $ dev.resids:function (y, mu, wt)  
 $ aic       :function (y, n, mu, wt, dev)  
 $ mu.eta    :function (eta)  
 $ initialize:  expression({  if (any(y <= 0))  stop("non-positive values not allowed for the 'Gamma' family")  n <- rep.int(1, n| __truncated__
 $ validmu   :function (mu)  
 $ valideta  :function (eta)  
 $ simulate  :function (object, nsim)  
 - attr(*, "class")= chr "family"

gf$linkinv

function (eta) 
1/eta
<environment: namespace:stats>

gf$variance(-3:4) #- == (.)^2

[1]  9  4  1  0  1  4  9 16

## Binomial with default 'logit' link:  Check some properties visually:
bi <- binomial()
et <- seq(-10,10, by=1/8)
plot(et, bi$mu.eta(et), type="l")
## show that mu.eta() is derivative of linkinv() :
lines((et[-1]+et[-length(et)])/2, col=adjustcolor("red", 1/4),
      diff(bi$linkinv(et))/diff(et), type="l", lwd=4)
## which here is the logistic density:
lines(et, dlogis(et), lwd=3, col=adjustcolor("blue", 1/4))
plot of chunk example-stats-family-1
stopifnot(exprs = {
  all.equal(bi$ mu.eta(et), dlogis(et))
  all.equal(bi$linkinv(et), plogis(et) -> m)
  all.equal(bi$linkfun(m ), qlogis(m))    #  logit(.) == qlogis(.) !
})

## Data from example(glm) :
d.AD <- data.frame(treatment = gl(3,3),
                   outcome   = gl(3,1,9),
                   counts    = c(18,17,15, 20,10,20, 25,13,12))
glm.D93 <- glm(counts ~ outcome + treatment, d.AD, family = poisson())
## Quasipoisson: compare with above / example(glm) :
glm.qD93 <- glm(counts ~ outcome + treatment, d.AD, family = quasipoisson())
## No test: 
glm.qD93

Call:  glm(formula = counts ~ outcome + treatment, family = quasipoisson(), 
    data = d.AD)

Coefficients:
(Intercept)     outcome2     outcome3   treatment2   treatment3  
  3.045e+00   -4.543e-01   -2.930e-01    1.338e-15    1.421e-15  

Degrees of Freedom: 8 Total (i.e. Null);  4 Residual
Null Deviance:	    10.58 
Residual Deviance: 5.129 	AIC: NA

anova  (glm.qD93, test = "F")

Analysis of Deviance Table

Model: quasipoisson, link: log

Response: counts

Terms added sequentially (first to last)


          Df Deviance Resid. Df Resid. Dev      F Pr(>F)
NULL                          8    10.5814              
outcome    2   5.4523         6     5.1291 2.1079  0.237
treatment  2   0.0000         4     5.1291 0.0000  1.000

summary(glm.qD93)

Call:
glm(formula = counts ~ outcome + treatment, family = quasipoisson(), 
    data = d.AD)

Deviance Residuals: 
       1         2         3         4         5         6         7         8  
-0.67125   0.96272  -0.16965  -0.21999  -0.95552   1.04939   0.84715  -0.09167  
       9  
-0.96656  

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.045e+00  1.944e-01  15.665  9.7e-05 ***
outcome2    -4.543e-01  2.299e-01  -1.976    0.119    
outcome3    -2.930e-01  2.192e-01  -1.337    0.252    
treatment2   1.338e-15  2.274e-01   0.000    1.000    
treatment3   1.421e-15  2.274e-01   0.000    1.000    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for quasipoisson family taken to be 1.2933)

    Null deviance: 10.5814  on 8  degrees of freedom
Residual deviance:  5.1291  on 4  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 4

## for Poisson results (same as from 'glm.D93' !) use
anova  (glm.qD93, dispersion = 1, test = "Chisq")

Analysis of Deviance Table

Model: quasipoisson, link: log

Response: counts

Terms added sequentially (first to last)


          Df Deviance Resid. Df Resid. Dev Pr(>Chi)  
NULL                          8    10.5814           
outcome    2   5.4523         6     5.1291  0.06547 .
treatment  2   0.0000         4     5.1291  1.00000  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(glm.qD93, dispersion = 1)

Call:
glm(formula = counts ~ outcome + treatment, family = quasipoisson(), 
    data = d.AD)

Deviance Residuals: 
       1         2         3         4         5         6         7         8  
-0.67125   0.96272  -0.16965  -0.21999  -0.95552   1.04939   0.84715  -0.09167  
       9  
-0.96656  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)  3.045e+00  1.709e-01  17.815   <2e-16 ***
outcome2    -4.543e-01  2.022e-01  -2.247   0.0246 *  
outcome3    -2.930e-01  1.927e-01  -1.520   0.1285    
treatment2   1.338e-15  2.000e-01   0.000   1.0000    
treatment3   1.421e-15  2.000e-01   0.000   1.0000    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for quasipoisson family taken to be 1)

    Null deviance: 10.5814  on 8  degrees of freedom
Residual deviance:  5.1291  on 4  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 4

## End(No test)


## Example of user-specified link, a logit model for p^days
## See Shaffer, T.  2004. Auk 121(2): 526-540.
logexp <- function(days = 1)
{
    linkfun <- function(mu) qlogis(mu^(1/days))
    linkinv <- function(eta) plogis(eta)^days
    mu.eta  <- function(eta) days * plogis(eta)^(days-1) *
                  binomial()$mu.eta(eta)
    valideta <- function(eta) TRUE
    link <- paste0("logexp(", days, ")")
    structure(list(linkfun = linkfun, linkinv = linkinv,
                   mu.eta = mu.eta, valideta = valideta, name = link),
              class = "link-glm")
}
(bil3 <- binomial(logexp(3)))

Family: binomial 
Link function: logexp(3) 

## Don't show: 
stopifnot(length(bil3$mu.eta(as.double(0:5))) == 6)
## End(Don't show)
## in practice this would be used with a vector of 'days', in
## which case use an offset of 0 in the corresponding formula
## to get the null deviance right.

## Binomial with identity link: often not a good idea, as both
## computationally and conceptually difficult:
binomial(link = "identity")  ## is exactly the same as

Family: binomial 
Link function: identity 

binomial(link = make.link("identity"))

Family: binomial 
Link function: identity 

## tests of quasi
x <- rnorm(100)
y <- rpois(100, exp(1+x))
glm(y ~ x, family = quasi(variance = "mu", link = "log"))

Call:  glm(formula = y ~ x, family = quasi(variance = "mu", link = "log"))

Coefficients:
(Intercept)            x  
     0.9988       0.9829  

Degrees of Freedom: 99 Total (i.e. Null);  98 Residual
Null Deviance:	    534.7 
Residual Deviance: 110.5 	AIC: NA

# which is the same as
glm(y ~ x, family = poisson)

Call:  glm(formula = y ~ x, family = poisson)

Coefficients:
(Intercept)            x  
     0.9988       0.9829  

Degrees of Freedom: 99 Total (i.e. Null);  98 Residual
Null Deviance:	    534.7 
Residual Deviance: 110.5 	AIC: 375.5

glm(y ~ x, family = quasi(variance = "mu^2", link = "log"))

Call:  glm(formula = y ~ x, family = quasi(variance = "mu^2", link = "log"))

Coefficients:
(Intercept)            x  
     1.0108       0.9649  

Degrees of Freedom: 99 Total (i.e. Null);  98 Residual
Null Deviance:	    86.3 
Residual Deviance: 24.2 	AIC: NA

## Not run: glm(y ~ x, family = quasi(variance = "mu^3", link = "log")) # fails
y <- rbinom(100, 1, plogis(x))
# need to set a starting value for the next fit
glm(y ~ x, family = quasi(variance = "mu(1-mu)", link = "logit"), start = c(0,1))

Call:  glm(formula = y ~ x, family = quasi(variance = "mu(1-mu)", link = "logit"), 
    start = c(0, 1))

Coefficients:
(Intercept)            x  
    -0.1826       0.7715  

Degrees of Freedom: 99 Total (i.e. Null);  98 Residual
Null Deviance:	    138 
Residual Deviance: 124.3 	AIC: NA
[Package stats version 4.2.3 Index]