Fazendo o gráfico.
x1 = c(1250,1300,1350,1250,1300,1250,1300,1350,1350)
x2 = c(6,7,6,7,6,8,8,7,8)
y = c(80, 95, 101, 85, 92, 87, 96, 106, 108)
library(scatterplot3d)
gr <- scatterplot3d(x1,
x2,
y)
gr
## $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: 0x505d4e0>
## <environment: 0x69d7238>
##
## $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: 0x31fdcd0>
## <environment: 0x69d7238>
##
## $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: 0x419a5e8>
## <environment: 0x69d7238>
##
## $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: 0x696e0e8>
## <environment: 0x69d7238>
##
## $contour3d
## function (f, x.count = 10, y.count = 10, type = "l", lty = "24",
## x.resolution = 50, y.resolution = 50, ...)
## {
## if (class(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 (class(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 (class(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: 0x6ac9a68>
## <environment: 0x69d7238>
##
## $par.mar
## $par.mar$mar
## [1] 5.1 4.1 4.1 2.1
De acordo com o gráfico, um modelo do primeiro grau seria o suficiente para ajustar um modelo aos dados observados.
Vamos ao ajuste do primeiro grau com a interação entre x1 e x2.
reg = lm(y ~ x1 + x2 + x1:x2)
summary(reg)
##
## Call:
## lm(formula = y ~ x1 + x2 + x1:x2)
##
## Residuals:
## 1 2 3 4 5 6 7 8
## -0.94444 0.55556 -0.94444 1.05556 0.55556 0.05556 -1.44444 1.05556
## 9
## 0.05556
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.996e+02 1.063e+02 -1.878 0.119
## x1 2.100e-01 8.172e-02 2.570 0.050 .
## x2 3.000e+00 1.508e+01 0.199 0.850
## x1:x2 -1.199e-16 1.160e-02 0.000 1.000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.16 on 5 degrees of freedom
## Multiple R-squared: 0.9907, Adjusted R-squared: 0.9851
## F-statistic: 177.4 on 3 and 5 DF, p-value: 1.697e-05
Percebe-se acima, que a interação entre x1 e x2 é não significativa. Logo,podemos retirar tal termo do modelo. Segue novo ajste.
reg1 = lm(y ~ x1 + x2)
summary(reg1)
##
## Call:
## lm(formula = y ~ x1 + x2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.44444 -0.94444 0.05556 0.55556 1.05556
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.996e+02 1.164e+01 -17.143 2.52e-06 ***
## x1 2.100e-01 8.642e-03 24.299 3.19e-07 ***
## x2 3.000e+00 4.321e-01 6.943 0.000443 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.058 on 6 degrees of freedom
## Multiple R-squared: 0.9907, Adjusted R-squared: 0.9876
## F-statistic: 319.3 on 2 and 6 DF, p-value: 8.064e-07
Percebemos nos resultados acima, que todos os coeficientes foram significativos (P<0,05) e de acordo com o coeficiente de determinação (\(R^2 = 0,9876\)) explica muito bem a variabilidade dos dados.
Apresentando um gráfico estático primeiro.
gr <- scatterplot3d(x1,
x2,
y)
gr$plane3d(reg1)
Apresentando um gráfico dinâmico.
# Fazendo o gráfico dinamicamente!
library(rgl)
y_estimado <- outer(X = x1, Y = x2, FUN = function(x1, x2) {
predict(reg1, newdata = data.frame(x1 = x1, x2 = x2))
})
dat <- data.frame(y,x1,x2)
with(data = dat,
plot3d(x = x1,
y = x2,
z = y_estimado,
type = "s", # s_phere
col = "red",
size = 2
)
)
## Adicionando a superfície ajustada
surface3d(x = x1,
y = x2,
z = y_estimado,
col = 3, alpha = 0.5)
