R version 4.5.0 beta (2025-03-27 r88065 ucrt) -- "How About a Twenty-Six" Copyright (C) 2025 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > require(glarma) Loading required package: glarma > ### Extract model from a glarma fit to use in simulation > ### Example with asthma data > data(Asthma) > y <- Asthma[,1] > X <- as.matrix(Asthma[,2:16]) > ## Pearson Residuals, Fisher Scoring > glarmamod <- glarma(y, X, thetaLags = 7, type = "Poi", method = "FS", + residuals = "Pearson", maxit = 100, grad = 1e-6) > mod <- extractGlarmaSimModel(glarmamod) > str(mod) List of 8 $ X : num [1:1461, 1:15] 1 1 1 1 1 1 1 1 1 1 ... ..- attr(*, "dimnames")=List of 2 .. ..$ : NULL .. ..$ : chr [1:15] "Intercept" "Sunday" "Monday" "CosAnnual" ... $ beta : Named num [1:15] 0.583 0.197 0.23 -0.214 0.176 ... ..- attr(*, "names")= chr [1:15] "Intercept" "Sunday" "Monday" "CosAnnual" ... $ type : chr "Poi" $ residType: chr "Pearson" $ phiLags : num(0) $ phi : num(0) $ thetaLags: num 7 $ theta : Named num 0.0423 ..- attr(*, "names")= chr "theta_7" - attr(*, "class")= chr "glarmaSimModel" > > glarmamod <- glarma(y, X, thetaLags = 7, type = "NegBin", method = "NR", + residuals = "Pearson", alphaInit = 0, + maxit = 100, grad = 1e-6) > mod <- extractGlarmaSimModel(glarmamod) > str(mod) List of 9 $ X : num [1:1461, 1:15] 1 1 1 1 1 1 1 1 1 1 ... ..- attr(*, "dimnames")=List of 2 .. ..$ : NULL .. ..$ : chr [1:15] "Intercept" "Sunday" "Monday" "CosAnnual" ... $ beta : Named num [1:15] 0.584 0.195 0.23 -0.215 0.177 ... ..- attr(*, "names")= chr [1:15] "Intercept" "Sunday" "Monday" "CosAnnual" ... $ type : chr "NegBin" $ alpha : Named num 37.2 ..- attr(*, "names")= chr "alpha" $ residType: chr "Pearson" $ phiLags : num(0) $ phi : num(0) $ thetaLags: num 7 $ theta : Named num 0.0439 ..- attr(*, "names")= chr "theta_7" - attr(*, "class")= chr "glarmaSimModel" > > data(DriverDeaths) > y <- DriverDeaths[, "Deaths"] > X <- as.matrix(DriverDeaths[, 2:5]) > Population <- DriverDeaths[, "Population"] > > ### Offset included > glarmamodOffset <- glarma(y, X, offset = log(Population/100000), + phiLags = c(12), + type = "Poi", method = "FS", + residuals = "Pearson", maxit = 100, grad = 1e-6) > mod <- extractGlarmaSimModel(glarmamodOffset) > str(mod) List of 9 $ X : num [1:72, 1:4] 1 1 1 1 1 1 1 1 1 1 ... ..- attr(*, "dimnames")=List of 2 .. ..$ : NULL .. ..$ : chr [1:4] "Intercept" "ReducedBAC" "FriSat" "lnOMVDRate" $ beta : Named num [1:4] -2.1352 -0.3859 0.0866 0.5572 ..- attr(*, "names")= chr [1:4] "Intercept" "ReducedBAC" "FriSat" "lnOMVDRate" $ offset : num [1:72] 2.13 2.13 2.13 2.13 2.13 ... $ type : chr "Poi" $ residType: chr "Pearson" $ phiLags : num 12 $ phi : Named num -0.0338 ..- attr(*, "names")= chr "phi_12" $ thetaLags: num(0) $ theta : num(0) - attr(*, "class")= chr "glarmaSimModel" > > data(OxBoatRace) > > y1 <- OxBoatRace$Camwin > n1 <- rep(1, length(OxBoatRace$Year)) > Y <- cbind(y1, n1 - y1) > X <- cbind(OxBoatRace$Intercept, OxBoatRace$Diff) > colnames(X) <- c("Intercept", "Weight Diff") > > oxcamglm <- glm(Y ~ Diff + I(Diff^2), + data = OxBoatRace, + family = binomial(link = "logit"), x = TRUE) > summary(oxcamglm) Call: glm(formula = Y ~ Diff + I(Diff^2), family = binomial(link = "logit"), data = OxBoatRace, x = TRUE) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.382472 0.202334 1.89 0.05872 . Diff 0.114623 0.037331 3.07 0.00214 ** I(Diff^2) -0.010097 0.004855 -2.08 0.03756 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 216.16 on 155 degrees of freedom Residual deviance: 198.62 on 153 degrees of freedom AIC: 204.62 Number of Fisher Scoring iterations: 5 > > X <- oxcamglm$x > > glarmamod <- glarma(Y, X, thetaLags = c(1, 2), type = "Bin", method = "NR", + residuals = "Pearson", maxit = 100, grad = 1e-6) > mod <- extractGlarmaSimModel(glarmamod) > str(mod) List of 9 $ X : num [1:156, 1:3] 1 1 1 1 1 1 1 1 1 1 ... ..- attr(*, "dimnames")=List of 2 .. ..$ : chr [1:156] "1" "2" "3" "4" ... .. ..$ : chr [1:3] "(Intercept)" "Diff" "I(Diff^2)" ..- attr(*, "assign")= int [1:3] 0 1 2 $ beta : Named num [1:3] 0.35 0.1148 -0.0113 ..- attr(*, "names")= chr [1:3] "(Intercept)" "Diff" "I(Diff^2)" $ type : chr "Bin" $ mBin : num [1:156] 1 1 1 1 1 1 1 1 1 1 ... $ residType: chr "Pearson" $ phiLags : num(0) $ phi : num(0) $ thetaLags: num [1:2] 1 2 $ theta : Named num [1:2] 0.34 0.555 ..- attr(*, "names")= chr [1:2] "theta_1" "theta_2" - attr(*, "class")= chr "glarmaSimModel" > > > data(RobberyConvict) > datalen <- dim(RobberyConvict)[1] > monthmat <- matrix(0, nrow = datalen, ncol = 12) > dimnames(monthmat) <- list(NULL, c("Jan", "Feb", "Mar", "Apr", "May", "Jun", + "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")) > months <- unique(months(strptime(RobberyConvict$Date, format = "%m/%d/%Y"), + abbreviate=TRUE)) > > > for (j in 1:12) { + monthmat[months(strptime(RobberyConvict$Date, "%m/%d/%Y"), + abbreviate = TRUE) == months[j], j] <- 1 + } > > RobberyConvict <- cbind(rep(1, datalen), RobberyConvict, monthmat) > rm(monthmat) > > ## LOWER COURT ROBBERY > y1 <- RobberyConvict$LC.Y > n1 <- RobberyConvict$LC.N > > Y <- cbind(y1, n1-y1) > > glm.LCRobbery <- glm(Y ~ Step.2001 + + I(Feb + Mar + Apr + May + Jun + Jul) + + I(Aug + Sep + Oct + Nov + Dec), + data = RobberyConvict, family = binomial(link = logit), + na.action = na.omit, x = TRUE) > > summary(glm.LCRobbery, corr = FALSE) Call: glm(formula = Y ~ Step.2001 + I(Feb + Mar + Apr + May + Jun + Jul) + I(Aug + Sep + Oct + Nov + Dec), family = binomial(link = logit), data = RobberyConvict, na.action = na.omit, x = TRUE) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.25685 0.15605 -1.646 0.09978 . Step.2001 0.82315 0.08135 10.119 < 2e-16 *** I(Feb + Mar + Apr + May + Jun + Jul) -0.37228 0.16188 -2.300 0.02146 * I(Aug + Sep + Oct + Nov + Dec) -0.50068 0.16549 -3.025 0.00248 ** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 327.48 on 149 degrees of freedom Residual deviance: 212.12 on 146 degrees of freedom AIC: 684.79 Number of Fisher Scoring iterations: 4 > > X <- glm.LCRobbery$x > > > ## Newton Raphson > glarmamod <- glarma(Y, X, phiLags = c(1), type = "Bin", method = "NR", + residuals = "Pearson", maxit = 100, grad = 1e-6) > mod <- extractGlarmaSimModel(glarmamod) > str(mod) List of 9 $ X : num [1:150, 1:4] 1 1 1 1 1 1 1 1 1 1 ... ..- attr(*, "dimnames")=List of 2 .. ..$ : chr [1:150] "1" "2" "3" "4" ... .. ..$ : chr [1:4] "(Intercept)" "Step.2001" "I(Feb + Mar + Apr + May + Jun + Jul)" "I(Aug + Sep + Oct + Nov + Dec)" ..- attr(*, "assign")= int [1:4] 0 1 2 3 $ beta : Named num [1:4] -0.275 0.822 -0.357 -0.5 ..- attr(*, "names")= chr [1:4] "(Intercept)" "Step.2001" "I(Feb + Mar + Apr + May + Jun + Jul)" "I(Aug + Sep + Oct + Nov + Dec)" $ type : chr "Bin" $ mBin : int [1:150] 12 11 15 15 11 17 12 15 16 10 ... $ residType: chr "Pearson" $ phiLags : num 1 $ phi : Named num 0.0818 ..- attr(*, "names")= chr "phi_1" $ thetaLags: num(0) $ theta : num(0) - attr(*, "class")= chr "glarmaSimModel" > > proc.time() user system elapsed 4.25 0.50 4.75