gnm-package: Generalized Nonlinear Models

Description Details Author(s) References See Also Examples

Description

Functions to specify, fit and evaluate generalized nonlinear models.

Details

gnm provides functions to fit generalized nonlinear models by maximum likelihood. Such models extend the class of generalized linear models by allowing nonlinear terms in the predictor.

Some special cases are models with multiplicative interaction terms, such as the UNIDIFF and row-column association models from sociology and the AMMI and GAMMI models from crop science; stereotype models for ordered categorical response, and diagonal reference models for dependence on a square two-way classification.

gnm is a major re-working of an earlier Xlisp-Stat package, "Llama". Over-parameterized representations of models are used throughout; functions are provided for inference on estimable parameter combinations, as well as standard methods for diagnostics etc.

The following documentation provides further information on the gnm package:

gnmOverview

vignette("gnmOverview", package = "gnm")

NEWS

file.show(system.file("NEWS", package = "gnm"))

Author(s)

Heather Turner and David Firth

Maintainer: Heather Turner <[email protected]>

References

http://www.warwick.ac.uk/go/gnm

See Also

gnm for the model fitting function, with links to associated functions.

Examples

1

Example output

	demo(gnm)
	---- ~~~

> message("1. Set seed as gnm returns random parameterization")
1. Set seed as gnm returns random parameterization

> set.seed(1)

> {
+     if (interactive()) {
+         cat("\n3. Type <Return> to fit (linear) uniform association model,  ",
+             "\n   using Diag() to fit diagonal effects: ")
+         readline()
+     }
+     else
+         message("2. Fit (linear) uniform association model, using Diag() to fit",
+                 "   diagonal effects")
+ }
2. Fit (linear) uniform association model, using Diag() to fit   diagonal effects

> Rscore <- scale(as.numeric(row(occupationalStatus)), scale = FALSE)

> Cscore <- scale(as.numeric(col(occupationalStatus)), scale = FALSE)

> Uniform <- gnm(Freq ~ origin + destination + Diag(origin, destination) +
+                Rscore:Cscore, family = poisson, data = occupationalStatus)

> summary(Uniform)

Call:
gnm(formula = Freq ~ origin + destination + Diag(origin, destination) + 
    Rscore:Cscore, family = poisson, data = occupationalStatus)

Deviance Residuals: 
       Min          1Q      Median          3Q         Max  
-2.652e+00  -6.267e-01  -1.054e-08   5.913e-01   2.096e+00  

Coefficients:
                           Estimate Std. Error z value Pr(>|z|)    
(Intercept)                0.568592   0.183358   3.101 0.001929 ** 
origin2                    0.431314   0.149415   2.887 0.003893 ** 
origin3                    1.461862   0.131141  11.147  < 2e-16 ***
origin4                    1.788731   0.126588  14.130  < 2e-16 ***
origin5                    0.441011   0.144844   3.045 0.002329 ** 
origin6                    2.491735   0.121219  20.556  < 2e-16 ***
origin7                    1.127536   0.129032   8.738  < 2e-16 ***
origin8                    0.796445   0.131863   6.040 1.54e-09 ***
destination2               0.873202   0.166902   5.232 1.68e-07 ***
destination3               1.813992   0.153861  11.790  < 2e-16 ***
destination4               2.082515   0.150887  13.802  < 2e-16 ***
destination5               1.366383   0.155590   8.782  < 2e-16 ***
destination6               2.816369   0.146054  19.283  < 2e-16 ***
destination7               1.903918   0.147810  12.881  < 2e-16 ***
destination8               1.398585   0.151658   9.222  < 2e-16 ***
Diag(origin, destination)1 1.665495   0.237383   7.016 2.28e-12 ***
Diag(origin, destination)2 0.959681   0.212122   4.524 6.06e-06 ***
Diag(origin, destination)3 0.021750   0.156530   0.139 0.889490    
Diag(origin, destination)4 0.226399   0.124364   1.820 0.068689 .  
Diag(origin, destination)5 0.808646   0.229754   3.520 0.000432 ***
Diag(origin, destination)6 0.132277   0.077191   1.714 0.086597 .  
Diag(origin, destination)7 0.506709   0.115936   4.371 1.24e-05 ***
Diag(origin, destination)8 0.221880   0.134803   1.646 0.099771 .  
Rscore:Cscore              0.136974   0.007489  18.289  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

Residual deviance: 58.436 on 40 degrees of freedom
AIC: 428.78

Number of iterations: 4


> {
+     if (interactive()) {
+         cat("\n3. Type <Return> to fit an association model using Mult() to fit",
+             "\n   separate row and column effects:")
+         readline()
+     }
+     else message("3. Fit an association model using Mult() to fit separate row and ",
+                  "column effects")
+ }
3. Fit an association model using Mult() to fit separate row and column effects

> RC <- gnm(Freq ~ origin + destination + Diag(origin, destination) +
+           Mult(origin, destination), family = poisson,
+           data = occupationalStatus)
Initialising
Running start-up iterations..
Running main iterations.........
Done

> summary(RC)

Call:
gnm(formula = Freq ~ origin + destination + Diag(origin, destination) + 
    Mult(origin, destination), family = poisson, data = occupationalStatus)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.6436  -0.3963   0.0000   0.3785   1.4915  

Coefficients:
                             Estimate Std. Error z value Pr(>|z|)   
(Intercept)                   0.11881         NA      NA       NA   
origin2                       0.49005         NA      NA       NA   
origin3                       1.74544         NA      NA       NA   
origin4                       2.15853         NA      NA       NA   
origin5                       0.93562         NA      NA       NA   
origin6                       3.09635         NA      NA       NA   
origin7                       1.83005         NA      NA       NA   
origin8                       1.60214         NA      NA       NA   
destination2                  0.95334         NA      NA       NA   
destination3                  1.87412         NA      NA       NA   
destination4                  2.04438         NA      NA       NA   
destination5                  1.45485         NA      NA       NA   
destination6                  2.81298         NA      NA       NA   
destination7                  1.88522         NA      NA       NA   
destination8                  1.42813         NA      NA       NA   
Diag(origin, destination)1    1.47923    0.45401   3.258  0.00112 **
Diag(origin, destination)2    0.47718    0.34175   1.396  0.16263   
Diag(origin, destination)3   -0.01813    0.18113  -0.100  0.92028   
Diag(origin, destination)4    0.38985    0.12743   3.059  0.00222 **
Diag(origin, destination)5    0.77988    0.23722   3.288  0.00101 **
Diag(origin, destination)6    0.13316    0.07922   1.681  0.09280 . 
Diag(origin, destination)7    0.44256    0.15348   2.884  0.00393 **
Diag(origin, destination)8    0.40731    0.21930   1.857  0.06327 . 
Mult(., destination).origin1 -1.74022         NA      NA       NA   
Mult(., destination).origin2 -1.56580         NA      NA       NA   
Mult(., destination).origin3 -0.67881         NA      NA       NA   
Mult(., destination).origin4  0.05605         NA      NA       NA   
Mult(., destination).origin5  0.27485         NA      NA       NA   
Mult(., destination).origin6  0.72410         NA      NA       NA   
Mult(., destination).origin7  1.34206         NA      NA       NA   
Mult(., destination).origin8  1.65900         NA      NA       NA   
Mult(origin, .).destination1 -1.32971         NA      NA       NA   
Mult(origin, .).destination2 -1.05346         NA      NA       NA   
Mult(origin, .).destination3 -0.66904         NA      NA       NA   
Mult(origin, .).destination4 -0.19797         NA      NA       NA   
Mult(origin, .).destination5 -0.25574         NA      NA       NA   
Mult(origin, .).destination6  0.21526         NA      NA       NA   
Mult(origin, .).destination7  0.51130         NA      NA       NA   
Mult(origin, .).destination8  0.66730         NA      NA       NA   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

Std. Error is NA where coefficient has been constrained or is unidentified

Residual deviance: 29.149 on 28 degrees of freedom
AIC: 423.49

Number of iterations: 9


> {
+     if (interactive()) {
+         cat("\n4. Type <Return> to fit an association model using MultHomog()",
+             "\n   to fit homogeneous row-column effects:")
+         readline()
+     }
+     else message("4. Fit an association model using MultHomog()\n",
+                  "to fit homogeneous row-column effects")
+ }
4. Fit an association model using MultHomog()
to fit homogeneous row-column effects

> RChomog <- gnm(Freq ~ origin + destination + Diag(origin, destination) +
+                MultHomog(origin, destination), family = poisson,
+                data = occupationalStatus)
Initialising
Running start-up iterations..
Running main iterations........
Done

> summary(RChomog)

Call:
gnm(formula = Freq ~ origin + destination + Diag(origin, destination) + 
    MultHomog(origin, destination), family = poisson, data = occupationalStatus)


Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.6588  -0.4297   0.0000   0.3862   1.7208  

Coefficients:
                                Estimate Std. Error z value Pr(>|z|)    
(Intercept)                      0.22383         NA      NA       NA    
origin2                          0.51137         NA      NA       NA    
origin3                          1.59736         NA      NA       NA    
origin4                          1.89707         NA      NA       NA    
origin5                          0.67717         NA      NA       NA    
origin6                          2.72219         NA      NA       NA    
origin7                          1.38190         NA      NA       NA    
origin8                          1.11206         NA      NA       NA    
destination2                     0.93038         NA      NA       NA    
destination3                     1.94178         NA      NA       NA    
destination4                     2.18551         NA      NA       NA    
destination5                     1.57658         NA      NA       NA    
destination6                     3.02567         NA      NA       NA    
destination7                     2.13350         NA      NA       NA    
destination8                     1.68743         NA      NA       NA    
Diag(origin, destination)1       1.52667    0.44658   3.419  0.00063 ***
Diag(origin, destination)2       0.45600    0.34595   1.318  0.18747    
Diag(origin, destination)3      -0.01598    0.18098  -0.088  0.92965    
Diag(origin, destination)4       0.38918    0.12748   3.053  0.00227 ** 
Diag(origin, destination)5       0.73852    0.23329   3.166  0.00155 ** 
Diag(origin, destination)6       0.13474    0.07934   1.698  0.08945 .  
Diag(origin, destination)7       0.45764    0.15103   3.030  0.00245 ** 
Diag(origin, destination)8       0.38847    0.22172   1.752  0.07976 .  
MultHomog(origin, destination)1  1.47021         NA      NA       NA    
MultHomog(origin, destination)2  1.25192         NA      NA       NA    
MultHomog(origin, destination)3  0.65375         NA      NA       NA    
MultHomog(origin, destination)4  0.06987         NA      NA       NA    
MultHomog(origin, destination)5  0.05270         NA      NA       NA    
MultHomog(origin, destination)6 -0.45905         NA      NA       NA    
MultHomog(origin, destination)7 -0.87520         NA      NA       NA    
MultHomog(origin, destination)8 -1.11877         NA      NA       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

Std. Error is NA where coefficient has been constrained or is unidentified

Residual deviance: 32.561 on 34 degrees of freedom
AIC: 414.9

Number of iterations: 8


> {
+     if (interactive()) {
+         cat("\n5. Type <Return> to compare models using anova:")
+         readline()
+     }
+     else message("5. Compare models using anova")
+ }
5. Compare models using anova

> anova(Uniform, RChomog, RC)
Analysis of Deviance Table

Model 1: Freq ~ origin + destination + Diag(origin, destination) + Rscore:Cscore
Model 2: Freq ~ origin + destination + Diag(origin, destination) + MultHomog(origin, 
    destination)
Model 3: Freq ~ origin + destination + Diag(origin, destination) + Mult(origin, 
    destination)
  Resid. Df Resid. Dev Df Deviance
1        40     58.436            
2        34     32.561  6  25.8748
3        28     29.149  6   3.4118

> message("6. Produce diagnostic plots for RChomog")
6. Produce diagnostic plots for RChomog

> plot(RChomog)

> message("7. Get simple constrasts of homogeneous row-column effects")
7. Get simple constrasts of homogeneous row-column effects

> getContrasts(RChomog, grep("MultHomog", names(coef(RChomog))))
                                  estimate        SE    quasiSE    quasiVar
MultHomog(origin, destination)1  0.0000000 0.0000000 0.15725493 0.024729113
MultHomog(origin, destination)2 -0.2182942 0.2346487 0.11901273 0.014164029
MultHomog(origin, destination)3 -0.8164601 0.1666550 0.06112070 0.003735739
MultHomog(origin, destination)4 -1.4003402 0.1602931 0.05183869 0.002687249
MultHomog(origin, destination)5 -1.4175083 0.1719512 0.07979360 0.006367019
MultHomog(origin, destination)6 -1.9292655 0.1573702 0.03598166 0.001294680
MultHomog(origin, destination)7 -2.3454075 0.1726343 0.07959733 0.006335734
MultHomog(origin, destination)8 -2.5889810 0.1886879 0.10953454 0.011997816

> message("End of demo. \n",
+  "See vignette(\"gnmOverview\", package = \"gnm\") for full manual.")
End of demo. 
See vignette("gnmOverview", package = "gnm") for full manual.
Warning messages:
1: In dropInf(r.w/(s * sqrt(1 - hii))) :
  Not plotting observations with leverage one:
  1, 10, 19, 28, 37, 46, 55, 64
2: In setCoefs - ref %*% setCoefs :
  Recycling array of length 1 in vector-array arithmetic is deprecated.
  Use c() or as.vector() instead.

gnm documentation built on June 24, 2018, 5 p.m.