Example 7.6.2.1 from Generalized Linear Mixed Models: Modern Concepts, Methods and Applications by Walter W. Stroup (p-231)
Source:R/Exam7.6.2.1.R
Exam7.6.2.1.RdExam7.6.2.1 Nonlinear Mean Models ( Quantitative by quantitative models)
References
Stroup, W. W. (2012). Generalized Linear Mixed Models: Modern Concepts, Methods and Applications. CRC Press.
@seealso
DataSet7.6
Author
Muhammad Yaseen (myaseen208@gmail.com)
Adeela Munawar (adeela.uaf@gmail.com)
Examples
library(scatterplot3d)
data(DataSet7.6)
library(dplyr)
#>
#> Attaching package: ‘dplyr’
#> The following object is masked from ‘package:MASS’:
#>
#> select
#> The following object is masked from ‘package:car’:
#>
#> recode
#> The following object is masked from ‘package:nlme’:
#>
#> collapse
#> The following objects are masked from ‘package:stats’:
#>
#> filter, lag
#> The following objects are masked from ‘package:base’:
#>
#> intersect, setdiff, setequal, union
library(magrittr)
DataSet7.6 <-
DataSet7.6 %>%
mutate(
logx1 = ifelse(test = x1 == 0, yes = log(x1 + 0.1), no = log(x1))
, logx2 = ifelse(test = x2 == 0, yes = log(x2 + 0.1), no = log(x2))
)
DataSet7.6
#> x1 x2 response logx1 logx2
#> 1 0 0 7.6 -2.3025851 -2.3025851
#> 2 0 1 6.9 -2.3025851 0.0000000
#> 3 0 2 8.1 -2.3025851 0.6931472
#> 4 0 3 8.7 -2.3025851 1.0986123
#> 5 0 4 7.7 -2.3025851 1.3862944
#> 6 1 0 7.4 0.0000000 -2.3025851
#> 7 1 1 14.1 0.0000000 0.0000000
#> 8 1 2 13.8 0.0000000 0.6931472
#> 9 1 3 13.8 0.0000000 1.0986123
#> 10 1 4 14.0 0.0000000 1.3862944
#> 11 2 0 7.6 0.6931472 -2.3025851
#> 12 2 1 13.3 0.6931472 0.0000000
#> 13 2 2 15.4 0.6931472 0.6931472
#> 14 2 3 14.8 0.6931472 1.0986123
#> 15 2 4 15.8 0.6931472 1.3862944
#> 16 3 0 6.1 1.0986123 -2.3025851
#> 17 3 1 11.7 1.0986123 0.0000000
#> 18 3 2 15.8 1.0986123 0.6931472
#> 19 3 3 17.1 1.0986123 1.0986123
#> 20 3 4 15.2 1.0986123 1.3862944
#> 21 4 0 5.7 1.3862944 -2.3025851
#> 22 4 1 14.4 1.3862944 0.0000000
#> 23 4 2 16.6 1.3862944 0.6931472
#> 24 4 3 17.0 1.3862944 1.0986123
#> 25 4 4 16.2 1.3862944 1.3862944
Exam7.6.2.1.lm <- lm(formula = response ~ x1*x2 + logx1*logx2 , data = DataSet7.6)
summary(Exam7.6.2.1.lm)
#>
#> Call:
#> lm(formula = response ~ x1 * x2 + logx1 * logx2, data = DataSet7.6)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -2.15804 -0.34365 0.05536 0.48754 1.45395
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 13.81950 1.04461 13.229 1.03e-10 ***
#> x1 -0.48366 0.44810 -1.079 0.294673
#> x2 -0.62134 0.44810 -1.387 0.182491
#> logx1 2.10830 0.33211 6.348 5.57e-06 ***
#> logx2 2.63409 0.33211 7.931 2.77e-07 ***
#> x1:x2 -0.06844 0.16211 -0.422 0.677871
#> logx1:logx2 0.84545 0.18510 4.568 0.000239 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
#> Residual standard error: 0.9269 on 18 degrees of freedom
#> Multiple R-squared: 0.9587, Adjusted R-squared: 0.945
#> F-statistic: 69.71 on 6 and 18 DF, p-value: 1.767e-11
#>
library(parameters)
model_parameters(Exam7.6.2.1.lm)
#> Parameter | Coefficient | SE | 95% CI | t(18) | p
#> --------------------------------------------------------------------
#> (Intercept) | 13.82 | 1.04 | [11.62, 16.01] | 13.23 | < .001
#> x1 | -0.48 | 0.45 | [-1.43, 0.46] | -1.08 | 0.295
#> x2 | -0.62 | 0.45 | [-1.56, 0.32] | -1.39 | 0.182
#> logx1 | 2.11 | 0.33 | [ 1.41, 2.81] | 6.35 | < .001
#> logx2 | 2.63 | 0.33 | [ 1.94, 3.33] | 7.93 | < .001
#> x1 × x2 | -0.07 | 0.16 | [-0.41, 0.27] | -0.42 | 0.678
#> logx1 × logx2 | 0.85 | 0.19 | [ 0.46, 1.23] | 4.57 | < .001
#>
#> Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed
#> using a Wald t-distribution approximation.
##---3D Scatter plot ( page#232)
attach(DataSet7.6)
(
ScatterPlot1 <-
scatterplot3d(
x = x1
, y = x2
, z = response
, color = response
, main = " 3D Scatter plot of response")
)
#> $xyz.convert
#> function (x, y = NULL, z = NULL)
#> {
#> xyz <- xyz.coords(x, y, z)
#> if (angle > 2) {
#> temp <- xyz$x
#> xyz$x <- xyz$y
#> xyz$y <- temp
#> }
#> y <- (xyz$y - y.add)/y.scal
#> return(list(x = xyz$x/x.scal + yx.f * y, y = xyz$z/z.scal +
#> yz.f * y))
#> }
#> <bytecode: 0x560506bf0b80>
#> <environment: 0x560506bda1d0>
#>
#> $points3d
#> function (x, y = NULL, z = NULL, type = "p", ...)
#> {
#> xyz <- xyz.coords(x, y, z)
#> if (angle > 2) {
#> temp <- xyz$x
#> xyz$x <- xyz$y
#> xyz$y <- temp
#> }
#> y2 <- (xyz$y - y.add)/y.scal
#> x <- xyz$x/x.scal + yx.f * y2
#> y <- xyz$z/z.scal + yz.f * y2
#> mem.par <- par(mar = mar, usr = usr)
#> if (type == "h") {
#> y2 <- z.min + yz.f * y2
#> segments(x, y, x, y2, ...)
#> points(x, y, type = "p", ...)
#> }
#> else points(x, y, type = type, ...)
#> }
#> <bytecode: 0x560506bee678>
#> <environment: 0x560506bda1d0>
#>
#> $plane3d
#> function (Intercept, x.coef = NULL, y.coef = NULL, lty = "dashed",
#> lty.box = NULL, draw_lines = TRUE, draw_polygon = FALSE,
#> polygon_args = list(border = NA, col = rgb(0, 0, 0, 0.2)),
#> ...)
#> {
#> if (!is.atomic(Intercept) && !is.null(coef(Intercept))) {
#> Intercept <- coef(Intercept)
#> if (!("(Intercept)" %in% names(Intercept)))
#> Intercept <- c(0, Intercept)
#> }
#> if (is.null(lty.box))
#> lty.box <- lty
#> if (is.null(x.coef) && length(Intercept) == 3) {
#> x.coef <- Intercept[if (angle > 2)
#> 3
#> else 2]
#> y.coef <- Intercept[if (angle > 2)
#> 2
#> else 3]
#> Intercept <- Intercept[1]
#> }
#> mem.par <- par(mar = mar, usr = usr)
#> x <- x.min:x.max
#> y <- 0:y.max
#> ltya <- c(lty.box, rep(lty, length(x) - 2), lty.box)
#> x.coef <- x.coef * x.scal
#> z1 <- (Intercept + x * x.coef + y.add * y.coef)/z.scal
#> z2 <- (Intercept + x * x.coef + (y.max * y.scal + y.add) *
#> y.coef)/z.scal
#> if (draw_polygon)
#> do.call("polygon", c(list(c(x.min, x.min + y.max * yx.f,
#> x.max + y.max * yx.f, x.max), c(z1[1], z2[1] + yz.f *
#> y.max, z2[length(z2)] + yz.f * y.max, z1[length(z1)])),
#> polygon_args))
#> if (draw_lines)
#> segments(x, z1, x + y.max * yx.f, z2 + yz.f * y.max,
#> lty = ltya, ...)
#> ltya <- c(lty.box, rep(lty, length(y) - 2), lty.box)
#> y.coef <- (y * y.scal + y.add) * y.coef
#> z1 <- (Intercept + x.min * x.coef + y.coef)/z.scal
#> z2 <- (Intercept + x.max * x.coef + y.coef)/z.scal
#> if (draw_lines)
#> segments(x.min + y * yx.f, z1 + y * yz.f, x.max + y *
#> yx.f, z2 + y * yz.f, lty = ltya, ...)
#> }
#> <bytecode: 0x560506bed3d0>
#> <environment: 0x560506bda1d0>
#>
#> $box3d
#> function (...)
#> {
#> mem.par <- par(mar = mar, usr = usr)
#> lines(c(x.min, x.max), c(z.max, z.max), ...)
#> lines(c(0, y.max * yx.f) + x.max, c(0, y.max * yz.f) + z.max,
#> ...)
#> lines(c(0, y.max * yx.f) + x.min, c(0, y.max * yz.f) + z.max,
#> ...)
#> lines(c(x.max, x.max), c(z.min, z.max), ...)
#> lines(c(x.min, x.min), c(z.min, z.max), ...)
#> lines(c(x.min, x.max), c(z.min, z.min), ...)
#> }
#> <bytecode: 0x560506be76e0>
#> <environment: 0x560506bda1d0>
#>
#> $contour3d
#> function (f, x.count = 10, y.count = 10, type = "l", lty = "24",
#> x.resolution = 50, y.resolution = 50, ...)
#> {
#> if (inherits(f, "lm")) {
#> vars <- all.vars(formula(f))
#> }
#> else vars <- c("z", "x", "y")
#> for (x1 in seq(x.range.fix[1], x.range.fix[2], length = x.count)) {
#> d <- data.frame(x1, seq(y.range.fix[1], y.range.fix[2],
#> length = y.resolution))
#> names(d) <- vars[-1]
#> if (inherits(f, "lm")) {
#> d[vars[1]] <- predict(f, newdata = d)
#> }
#> else d[vars[1]] <- f(d[[1]], d[[2]])
#> xyz <- xyz.coords(d)
#> if (angle > 2) {
#> temp <- xyz$x
#> xyz$x <- xyz$y
#> xyz$y <- temp
#> }
#> y2 <- (xyz$y - y.add)/y.scal
#> x <- xyz$x/x.scal + yx.f * y2
#> y <- xyz$z/z.scal + yz.f * y2
#> mem.par <- par(mar = mar, usr = usr)
#> if (type == "h") {
#> y2 <- z.min + yz.f * y2
#> segments(x, y, x, y2, ...)
#> points(x, y, type = "p", ...)
#> }
#> else points(x, y, type = type, lty = lty, ...)
#> }
#> for (x2 in seq(y.range.fix[1], y.range.fix[2], length = y.count)) {
#> d <- data.frame(seq(x.range.fix[1], x.range.fix[2], length = x.resolution),
#> x2)
#> names(d) <- vars[-1]
#> if (inherits(f, "lm")) {
#> d[vars[1]] <- predict(f, newdata = d)
#> }
#> else d[vars[1]] <- f(d[[1]], d[[2]])
#> xyz <- xyz.coords(d)
#> if (angle > 2) {
#> temp <- xyz$x
#> xyz$x <- xyz$y
#> xyz$y <- temp
#> }
#> y2 <- (xyz$y - y.add)/y.scal
#> x <- xyz$x/x.scal + yx.f * y2
#> y <- xyz$z/z.scal + yz.f * y2
#> mem.par <- par(mar = mar, usr = usr)
#> if (type == "h") {
#> y2 <- z.min + yz.f * y2
#> segments(x, y, x, y2, ...)
#> points(x, y, type = "p", ...)
#> }
#> else points(x, y, type = type, lty = lty, ...)
#> }
#> }
#> <bytecode: 0x560506be3d70>
#> <environment: 0x560506bda1d0>
#>
#> $par.mar
#> $par.mar$mar
#> [1] 5.1 4.1 4.1 2.1
#>
#>
##--- scatter plot with regression plane by using Hoerl function ( page#233)
grid.lines <- 5
x1.pred <- seq(min(x1), max(x1), length.out = grid.lines)
x2.pred <- seq(min(x2), max(x2), length.out = grid.lines)
x1x2 <- expand.grid( x = x1.pred, y = x2.pred)
z.pred <- matrix(data = predict(Exam7.6.2.1.lm, newdata = x1x2),
nrow = grid.lines
, ncol = grid.lines)
(ScatterPlot2 <-
scatterplot3d(
x = x1
, y = x2
, z = response
, pch = 20
, phi = 25
, theta = 30
, ticktype = "detailed"
, xlab = "x1"
, ylab = "x2"
, zlab = "response"
, add = FALSE
, surf = list(x = x1.pred ,
y = x2.pred ,
z = z.pred ,
facets = NA
)
, plot = TRUE
, main = "Fitted Response Surface by Hoerl Function"
)
)
#> Warning: "phi" is not a graphical parameter
#> Warning: "theta" is not a graphical parameter
#> Warning: "ticktype" is not a graphical parameter
#> Warning: "add" is not a graphical parameter
#> Warning: "surf" is not a graphical parameter
#> Warning: "plot" is not a graphical parameter
#> Warning: "phi" is not a graphical parameter
#> Warning: "theta" is not a graphical parameter
#> Warning: "ticktype" is not a graphical parameter
#> Warning: "add" is not a graphical parameter
#> Warning: "surf" is not a graphical parameter
#> Warning: "plot" is not a graphical parameter
#> $xyz.convert
#> function (x, y = NULL, z = NULL)
#> {
#> xyz <- xyz.coords(x, y, z)
#> if (angle > 2) {
#> temp <- xyz$x
#> xyz$x <- xyz$y
#> xyz$y <- temp
#> }
#> y <- (xyz$y - y.add)/y.scal
#> return(list(x = xyz$x/x.scal + yx.f * y, y = xyz$z/z.scal +
#> yz.f * y))
#> }
#> <bytecode: 0x560506bf0b80>
#> <environment: 0x56050459d9e0>
#>
#> $points3d
#> function (x, y = NULL, z = NULL, type = "p", ...)
#> {
#> xyz <- xyz.coords(x, y, z)
#> if (angle > 2) {
#> temp <- xyz$x
#> xyz$x <- xyz$y
#> xyz$y <- temp
#> }
#> y2 <- (xyz$y - y.add)/y.scal
#> x <- xyz$x/x.scal + yx.f * y2
#> y <- xyz$z/z.scal + yz.f * y2
#> mem.par <- par(mar = mar, usr = usr)
#> if (type == "h") {
#> y2 <- z.min + yz.f * y2
#> segments(x, y, x, y2, ...)
#> points(x, y, type = "p", ...)
#> }
#> else points(x, y, type = type, ...)
#> }
#> <bytecode: 0x560506bee678>
#> <environment: 0x56050459d9e0>
#>
#> $plane3d
#> function (Intercept, x.coef = NULL, y.coef = NULL, lty = "dashed",
#> lty.box = NULL, draw_lines = TRUE, draw_polygon = FALSE,
#> polygon_args = list(border = NA, col = rgb(0, 0, 0, 0.2)),
#> ...)
#> {
#> if (!is.atomic(Intercept) && !is.null(coef(Intercept))) {
#> Intercept <- coef(Intercept)
#> if (!("(Intercept)" %in% names(Intercept)))
#> Intercept <- c(0, Intercept)
#> }
#> if (is.null(lty.box))
#> lty.box <- lty
#> if (is.null(x.coef) && length(Intercept) == 3) {
#> x.coef <- Intercept[if (angle > 2)
#> 3
#> else 2]
#> y.coef <- Intercept[if (angle > 2)
#> 2
#> else 3]
#> Intercept <- Intercept[1]
#> }
#> mem.par <- par(mar = mar, usr = usr)
#> x <- x.min:x.max
#> y <- 0:y.max
#> ltya <- c(lty.box, rep(lty, length(x) - 2), lty.box)
#> x.coef <- x.coef * x.scal
#> z1 <- (Intercept + x * x.coef + y.add * y.coef)/z.scal
#> z2 <- (Intercept + x * x.coef + (y.max * y.scal + y.add) *
#> y.coef)/z.scal
#> if (draw_polygon)
#> do.call("polygon", c(list(c(x.min, x.min + y.max * yx.f,
#> x.max + y.max * yx.f, x.max), c(z1[1], z2[1] + yz.f *
#> y.max, z2[length(z2)] + yz.f * y.max, z1[length(z1)])),
#> polygon_args))
#> if (draw_lines)
#> segments(x, z1, x + y.max * yx.f, z2 + yz.f * y.max,
#> lty = ltya, ...)
#> ltya <- c(lty.box, rep(lty, length(y) - 2), lty.box)
#> y.coef <- (y * y.scal + y.add) * y.coef
#> z1 <- (Intercept + x.min * x.coef + y.coef)/z.scal
#> z2 <- (Intercept + x.max * x.coef + y.coef)/z.scal
#> if (draw_lines)
#> segments(x.min + y * yx.f, z1 + y * yz.f, x.max + y *
#> yx.f, z2 + y * yz.f, lty = ltya, ...)
#> }
#> <bytecode: 0x560506bed3d0>
#> <environment: 0x56050459d9e0>
#>
#> $box3d
#> function (...)
#> {
#> mem.par <- par(mar = mar, usr = usr)
#> lines(c(x.min, x.max), c(z.max, z.max), ...)
#> lines(c(0, y.max * yx.f) + x.max, c(0, y.max * yz.f) + z.max,
#> ...)
#> lines(c(0, y.max * yx.f) + x.min, c(0, y.max * yz.f) + z.max,
#> ...)
#> lines(c(x.max, x.max), c(z.min, z.max), ...)
#> lines(c(x.min, x.min), c(z.min, z.max), ...)
#> lines(c(x.min, x.max), c(z.min, z.min), ...)
#> }
#> <bytecode: 0x560506be76e0>
#> <environment: 0x56050459d9e0>
#>
#> $contour3d
#> function (f, x.count = 10, y.count = 10, type = "l", lty = "24",
#> x.resolution = 50, y.resolution = 50, ...)
#> {
#> if (inherits(f, "lm")) {
#> vars <- all.vars(formula(f))
#> }
#> else vars <- c("z", "x", "y")
#> for (x1 in seq(x.range.fix[1], x.range.fix[2], length = x.count)) {
#> d <- data.frame(x1, seq(y.range.fix[1], y.range.fix[2],
#> length = y.resolution))
#> names(d) <- vars[-1]
#> if (inherits(f, "lm")) {
#> d[vars[1]] <- predict(f, newdata = d)
#> }
#> else d[vars[1]] <- f(d[[1]], d[[2]])
#> xyz <- xyz.coords(d)
#> if (angle > 2) {
#> temp <- xyz$x
#> xyz$x <- xyz$y
#> xyz$y <- temp
#> }
#> y2 <- (xyz$y - y.add)/y.scal
#> x <- xyz$x/x.scal + yx.f * y2
#> y <- xyz$z/z.scal + yz.f * y2
#> mem.par <- par(mar = mar, usr = usr)
#> if (type == "h") {
#> y2 <- z.min + yz.f * y2
#> segments(x, y, x, y2, ...)
#> points(x, y, type = "p", ...)
#> }
#> else points(x, y, type = type, lty = lty, ...)
#> }
#> for (x2 in seq(y.range.fix[1], y.range.fix[2], length = y.count)) {
#> d <- data.frame(seq(x.range.fix[1], x.range.fix[2], length = x.resolution),
#> x2)
#> names(d) <- vars[-1]
#> if (inherits(f, "lm")) {
#> d[vars[1]] <- predict(f, newdata = d)
#> }
#> else d[vars[1]] <- f(d[[1]], d[[2]])
#> xyz <- xyz.coords(d)
#> if (angle > 2) {
#> temp <- xyz$x
#> xyz$x <- xyz$y
#> xyz$y <- temp
#> }
#> y2 <- (xyz$y - y.add)/y.scal
#> x <- xyz$x/x.scal + yx.f * y2
#> y <- xyz$z/z.scal + yz.f * y2
#> mem.par <- par(mar = mar, usr = usr)
#> if (type == "h") {
#> y2 <- z.min + yz.f * y2
#> segments(x, y, x, y2, ...)
#> points(x, y, type = "p", ...)
#> }
#> else points(x, y, type = type, lty = lty, ...)
#> }
#> }
#> <bytecode: 0x560506be3d70>
#> <environment: 0x56050459d9e0>
#>
#> $par.mar
#> $par.mar$mar
#> [1] 5.1 3.1 4.1 3.1
#>
#>