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. > ### Testing of glarmaSim with examples from package > require(glarma) Loading required package: glarma > ### Polio data > data("Polio") > y <- Polio[, 2] > X <- as.matrix(Polio[, 3:8]) > glarmamod <- glarma(y, X, thetaLags = c(1,2,5), type = "Poi", method = "FS", + residuals = "Pearson", maxit = 100, grad = 1e-6) > PolioModel <- extractGlarmaSimModel(glarmamod) > str(PolioModel) List of 8 $ X : num [1:168, 1:6] 1 1 1 1 1 1 1 1 1 1 ... ..- attr(*, "dimnames")=List of 2 .. ..$ : NULL .. ..$ : chr [1:6] "Intcpt" "Trend" "CosAnnual" "SinAnnual" ... $ beta : Named num [1:6] 0.13 -3.9284 -0.0991 -0.5308 0.2111 ... ..- attr(*, "names")= chr [1:6] "Intcpt" "Trend" "CosAnnual" "SinAnnual" ... $ type : chr "Poi" $ residType: chr "Pearson" $ phiLags : num(0) $ phi : num(0) $ thetaLags: num [1:3] 1 2 5 $ theta : Named num [1:3] 0.2185 0.1272 0.0873 ..- attr(*, "names")= chr [1:3] "theta_1" "theta_2" "theta_5" - attr(*, "class")= chr "glarmaSimModel" > par(mfrow = c(3,1)) > sim <- glarmaSim(PolioModel) > ts.plot(sim$W) > ts.plot(sim$mu) > ts.plot(sim$Y) > > ### Example with Oxford-Cambridge Boat Race > 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) > str(glarmamod) List of 27 $ delta : num [1:5, 1] 0.35 0.1148 -0.0113 0.3396 0.5552 ..- attr(*, "names")= chr [1:5] "(Intercept)" "Diff" "I(Diff^2)" "theta_1" ... $ logLik : num -91.5 $ logLikDeriv : num [1:5, 1] -3.89e-09 9.80e-09 -1.42e-07 3.25e-09 -2.17e-09 $ logLikDeriv2 : num [1:5, 1:5] -17.01 -1.48 -390.37 -2.02 -2 ... $ eta : num [1:156, 1] 0.34995 0.44169 -1.84625 -0.00776 0.49659 ... ..- attr(*, "dimnames")=List of 2 .. ..$ : chr [1:156] "1" "2" "3" "4" ... .. ..$ : NULL $ W : num [1:156(1d)] 0.35 0.0372 -2.1743 1.5442 2.3 ... $ residuals : num [1:156(1d)] -1.191 0.982 2.966 0.462 0.317 ... $ mu : num [1:156(1d)] 0.587 0.509 0.102 0.824 0.909 ... $ fitted.values: num [1:156(1d)] 0.587 0.509 0.102 0.824 0.909 ... $ cov : num [1:5, 1:5] 0.071422 0.001276 -0.000536 -0.001398 -0.002435 ... $ null.deviance: num 216 $ df.null : int 155 $ r : int 3 $ pq : int 2 $ phiLags : num(0) $ thetaLags : num [1:2] 1 2 $ y : num [1:156, 1:2] 0 1 1 1 1 0 1 1 1 0 ... ..- attr(*, "dimnames")=List of 2 .. ..$ : NULL .. ..$ : chr [1:2] "y1" "" $ 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 $ type : chr "Bin" $ method : chr "NR" $ residType : chr "Pearson" $ call : language glarma(y = Y, X = X, type = "Bin", method = "NR", residuals = "Pearson", thetaLags = c(1, 2), maxit = 100, grad = 1e-06) $ iter : num 5 $ errCode : num 0 $ WError : num 0 $ min : num 1.42e-07 $ aic : num 193 - attr(*, "class")= chr "glarma" > BoatRaceModel <- glarmaSimModel(X, beta = coef.glarma(glarmamod, type = "beta"), + phiLags = glarmamod$phiLags, + phi = rep(0.1, length(glarmamod$phiLags)), + thetaLags = glarmamod$thetaLags, + theta = c(0.34,0.56), + type = glarmamod$type, m = n1, + residType = glarmamod$residType) > str(BoatRaceModel) 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 : num [1:2] 0.34 0.56 - attr(*, "class")= chr "glarmaSimModel" > BoatRaceModel <- extractGlarmaSimModel(glarmamod) > str(BoatRaceModel) 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" > par(mfrow = c(3,1)) > sim <- glarmaSim(BoatRaceModel) > ts.plot(sim$W) > ts.plot(sim$mu) > ts.plot(sim$Y) > > > ### Example with asthma data, negative binomial > data(Asthma) > y <- Asthma[, 1] > X <- as.matrix(Asthma[, 2:16]) > > ## Pearson Residuals, Newton Raphson, Negative Binomial > ## Initial value of the shape parameter take to be zero > glarmamod <- glarma(y, X, thetaLags = 7, type = "NegBin", method = "NR", + residuals = "Pearson", alphaInit = 0, + maxit = 100, grad = 1e-6) > summary(glarmamod) Call: glarma(y = y, X = X, type = "NegBin", method = "NR", residuals = "Pearson", thetaLags = 7, alphaInit = 0, maxit = 100, grad = 1e-06) Pearson Residuals: Min 1Q Median 3Q Max -1.8491 -0.7406 -0.1754 0.6092 6.1776 Negative Binomial Parameter: Estimate Std.Error z-ratio Pr(>|z|) alpha 37.19 25.44 1.462 0.144 GLARMA Coefficients: Estimate Std.Error z-ratio Pr(>|z|) theta_7 0.04392 0.01936 2.269 0.0233 * Linear Model Coefficients: Estimate Std.Error z-ratio Pr(>|z|) Intercept 0.58397 0.06331 9.225 < 2e-16 *** Sunday 0.19455 0.05760 3.377 0.000732 *** Monday 0.22999 0.05642 4.076 4.57e-05 *** CosAnnual -0.21450 0.03965 -5.410 6.31e-08 *** SinAnnual 0.17728 0.04153 4.269 1.96e-05 *** H7 0.16843 0.05634 2.990 0.002794 ** NO2max -0.10404 0.03392 -3.067 0.002163 ** T1.1990 0.19903 0.05845 3.405 0.000662 *** T2.1990 0.13087 0.05897 2.219 0.026477 * T1.1991 0.08587 0.06746 1.273 0.203058 T2.1991 0.17082 0.05951 2.871 0.004096 ** T1.1992 0.25276 0.05669 4.459 8.24e-06 *** T2.1992 0.30572 0.05103 5.991 2.08e-09 *** T1.1993 0.43607 0.05233 8.334 < 2e-16 *** T2.1993 0.11412 0.06269 1.821 0.068679 . Null deviance: 1989.9 on 1460 degrees of freedom Residual deviance: 1442.6 on 1444 degrees of freedom AIC: 4873.511 Number of Newton Raphson iterations: 6 LRT and Wald Test: Alternative hypothesis: model is a GLARMA process Null hypothesis: model is a GLM with the same regression structure Statistic p-value LR Test 7.047 0.00794 ** Wald Test 5.147 0.02329 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > str(glarmamod) List of 27 $ delta : num [1:17, 1] 0.584 0.195 0.23 -0.215 0.177 ... ..- attr(*, "names")= chr [1:17] "Intercept" "Sunday" "Monday" "CosAnnual" ... $ logLik : num -2421 $ logLikDeriv : num [1:17, 1] -6.59e-09 -2.96e-08 7.77e-09 -7.17e-09 -5.83e-09 ... $ logLikDeriv2 : num [1:17, 1:17] -2368 -387 -393 153 -392 ... $ eta : num [1:1461, 1] 0.53 0.349 0.255 0.379 0.364 ... $ W : num [1:1461(1d)] 0.53 0.349 0.255 0.379 0.364 ... $ residuals : num [1:1461(1d)] 0.976 -0.345 0.614 -0.374 0.459 ... $ mu : num [1:1461(1d)] 1.7 1.42 1.29 1.46 1.44 ... $ fitted.values: num [1:1461(1d)] 1.7 1.42 1.29 1.46 1.44 ... $ cov : num [1:17, 1:17] 0.004008 -0.00118 -0.000775 -0.000897 0.000309 ... $ null.deviance: num 1990 $ df.null : int 1460 $ r : int 15 $ pq : int 1 $ phiLags : num(0) $ thetaLags : num 7 $ y : int [1:1461] 3 1 2 1 2 1 2 0 1 0 ... $ 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" ... $ type : chr "NegBin" $ method : chr "NR" $ residType : chr "Pearson" $ call : language glarma(y = y, X = X, type = "NegBin", method = "NR", residuals = "Pearson", thetaLags = 7, alphaInit = 0, ma| __truncated__ $ iter : num 6 $ errCode : num 0 $ WError : num 0 $ min : num 1.64e-07 $ aic : num 4874 - attr(*, "class")= chr "glarma" > AsthmaModel <- glarmaSimModel(X, beta = coef.glarma(glarmamod, type = "beta"), + phiLags = glarmamod$phiLags, + phi = rep(0.1, length(glarmamod$phiLags)), + thetaLags = glarmamod$thetaLags, + theta = c(0.044), + type = glarmamod$type, alpha = 37.19, + residType = glarmamod$residType) > str(AsthmaModel) 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 : num 37.2 $ residType: chr "Pearson" $ phiLags : num(0) $ phi : num(0) $ thetaLags: num 7 $ theta : num 0.044 - attr(*, "class")= chr "glarmaSimModel" > AsthmaModel <- extractGlarmaSimModel(glarmamod) > str(AsthmaModel) 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" > par(mfrow = c(3,1)) > sim <- glarmaSim(AsthmaModel) > ts.plot(sim$W) > ts.plot(sim$mu) > ts.plot(sim$Y) > > proc.time() user system elapsed 4.12 0.28 4.37