Rmd

rockchalk

Rmd

With library(rockchalk), graphs are not interactive but knittable.

doc

df <- read.delim("http://md.psych.bio.uni-goettingen.de/data/div/sat.txt")
m.i <- lm(univ_GPA ~ math_SAT * verb_SAT, data=df)

require(rockchalk)
## Loading required package: rockchalk
plotPlane(model = m.i, plotx1 = "math_SAT", plotx2 = "verb_SAT")

# or to emphasize the view of the interaction
plotPlane(model = m.i, plotx1 = "math_SAT", plotx2 = "verb_SAT", theta = 0, phi = 10)

plotPlane(model = m.i, plotx1 = "math_SAT", plotx2 = "verb_SAT", drawArrows = T)

Example Fishes

source

# create data frame df
df <- structure(list(Sample = structure(1:12, .Label = c("r.2002.c.1", 
  "r.2002.c.2", "r.2002.c.3", "r.2002.c.4", "r.2002.c.5", "r.2002.c.6", 
  "r.2002.s.1", "r.2002.s.2", "r.2002.s.3", "r.2002.s.4", "r.2002.s.5", 
  "r.2002.s.6"), class = "factor"), Coral.cover = c(59L, 71L, 55L, 
  58L, 62L, 65L, 86L, 71L, 91L, 63L, 75L, 63L), Coral.richness = c(5L, 
  9L, 10L, 10L, 8L, 10L, 10L, 10L, 6L, 0L, 2L, 2L), Water.Temp = c(23.9, 
  33.43, 23.33, 22.14, 26.82, 29.07, 38.76, 30.28, 27.85, 28.03, 
  31.67, 29.43), fish.rich = c(57L, 70L, 55L, 56L, 61L, 64L, 84L, 
  68L, 89L, 62L, 75L, 63L)), .Names = c("Sample", "Coral.cover", 
  "Coral.richness", "Water.Temp", "fish.rich"), class = "data.frame", row.names = c(NA, 
  -12L))

Coral.cover <- df[,2]
Coral.richness <- df[,3]
fish.rich <- df[,5]

# libraries
require(Hmisc)
## Loading required package: Hmisc
## Loading required package: lattice
## Loading required package: survival
## Loading required package: Formula
## Loading required package: ggplot2
## Registered S3 methods overwritten by 'ggplot2':
##   method         from 
##   [.quosures     rlang
##   c.quosures     rlang
##   print.quosures rlang
## 
## Attaching package: 'Hmisc'
## The following object is masked from 'package:rockchalk':
## 
##     summarize
## The following objects are masked from 'package:base':
## 
##     format.pval, units
require(car)
## Loading required package: car
## Loading required package: carData
## Registered S3 methods overwritten by 'car':
##   method                          from
##   influence.merMod                lme4
##   cooks.distance.influence.merMod lme4
##   dfbeta.influence.merMod         lme4
##   dfbetas.influence.merMod        lme4
# require(vegan)
require(rockchalk)
require(reshape)
## Loading required package: reshape
require(persp)
## Loading required package: persp
## Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
## logical.return = TRUE, : there is no package called 'persp'
mod3 <- mcGraph3(Coral.cover, Coral.richness, fish.rich, 
             interaction = F,
             theta = -40 ,
             phi =20,
             x1lab = "",
             x2lab = "",
             ylab = "",
             x1lim = c(46,100),
             ylim = c(-2.5, 12.5),
             zlim =c(40, 100),
             r=10,
             col = 'white',
             border = 'black',
             box = T,
             axes = T,
             ticktype = "detailed",
             ntick = 4)
## Warning in persp.default(x = c(46, 59.5, 73, 86.5, 100), y = c(-2.5,
## 1.25, : surface extends beyond the box

Version: 08 Juni, 2020 12:59

plot3d

Example Fishes

source

# create data frame df
df <- structure(list(Sample = structure(1:12, .Label = c("r.2002.c.1", 
  "r.2002.c.2", "r.2002.c.3", "r.2002.c.4", "r.2002.c.5", "r.2002.c.6", 
  "r.2002.s.1", "r.2002.s.2", "r.2002.s.3", "r.2002.s.4", "r.2002.s.5", 
  "r.2002.s.6"), class = "factor"), Coral.cover = c(59L, 71L, 55L, 
  58L, 62L, 65L, 86L, 71L, 91L, 63L, 75L, 63L), Coral.richness = c(5L, 
  9L, 10L, 10L, 8L, 10L, 10L, 10L, 6L, 0L, 2L, 2L), Water.Temp = c(23.9, 
  33.43, 23.33, 22.14, 26.82, 29.07, 38.76, 30.28, 27.85, 28.03, 
  31.67, 29.43), fish.rich = c(57L, 70L, 55L, 56L, 61L, 64L, 84L, 
  68L, 89L, 62L, 75L, 63L)), .Names = c("Sample", "Coral.cover", 
  "Coral.richness", "Water.Temp", "fish.rich"), class = "data.frame", row.names = c(NA, 
  -12L))

Coral.cover <- df[,2]
Coral.richness <- df[,3]
fish.rich <- df[,5]

## preparation of some values for mesh of fitted value
fit <- lm(fish.rich ~ Coral.cover + Coral.richness)  # model
x.p <- seq(46, 100, length = 20)                     # x-grid of mesh
y.p <- seq(-2.5, 12.5, length = 20)                  # y-grid of mesh
z.p <- matrix(predict(fit, expand.grid(Coral.cover = x.p, Coral.richness = y.p)), 20) # prediction from xy-grid

library(plot3D)
# box, grid, bottom points, and so on
scatter3D(Coral.cover, Coral.richness, rep(40, 12), colvar = NA, bty = "b2", 
          xlim = c(46,100), ylim = c(-2.5, 12.5), zlim = c(40,100), 
          # theta = -40, phi = 20, 
          theta =0, phi = 20, 
          r = 10, ticktype = "detailed", pch = 19, col = "gray", nticks = 4)
# mesh, real points
scatter3D(Coral.cover, Coral.richness, fish.rich, add = T, colvar = NA, col = "blue",
          surf = list(x = x.p, y = y.p, z = z.p, facets = NA, col = "gray80"))
# arrow from prediction to observation
arrows3D(x0 = Coral.cover, y0 = Coral.richness, z0 = fit$fitted.values, z1 = fish.rich, 
         type = "simple", lty = 2, add = T, col = "red")

base system persp

Example Fishes

source

# create data frame df
df <- structure(list(Sample = structure(1:12, .Label = c("r.2002.c.1", 
  "r.2002.c.2", "r.2002.c.3", "r.2002.c.4", "r.2002.c.5", "r.2002.c.6", 
  "r.2002.s.1", "r.2002.s.2", "r.2002.s.3", "r.2002.s.4", "r.2002.s.5", 
  "r.2002.s.6"), class = "factor"), Coral.cover = c(59L, 71L, 55L, 
  58L, 62L, 65L, 86L, 71L, 91L, 63L, 75L, 63L), Coral.richness = c(5L, 
  9L, 10L, 10L, 8L, 10L, 10L, 10L, 6L, 0L, 2L, 2L), Water.Temp = c(23.9, 
  33.43, 23.33, 22.14, 26.82, 29.07, 38.76, 30.28, 27.85, 28.03, 
  31.67, 29.43), fish.rich = c(57L, 70L, 55L, 56L, 61L, 64L, 84L, 
  68L, 89L, 62L, 75L, 63L)), .Names = c("Sample", "Coral.cover", 
  "Coral.richness", "Water.Temp", "fish.rich"), class = "data.frame", row.names = c(NA, 
  -12L))

Coral.cover <- df[,2]
Coral.richness <- df[,3]
fish.rich <- df[,5]

# libraries
# require(Hmisc)
# require(car)
# # require(vegan)
# require(rockchalk)
# require(reshape)
# require(persp)


## preparation of some values for mesh of fitted value
fit <- lm(fish.rich ~ Coral.cover + Coral.richness)  # model
x.p <- seq(46, 100, length = 20)                     # x-grid of mesh
y.p <- seq(-2.5, 12.5, length = 20)                  # y-grid of mesh
z.p <- matrix(predict(fit, expand.grid(Coral.cover = x.p, Coral.richness = y.p)), 20) # prediction from xy-grid


### [bonus] persp() version
pmat <- persp(x.p, y.p, z.p, xlim = c(46,100), ylim = c(-2.5, 12.5), zlim = c(40,100), theta = -40, phi = 20, 
              r = 10, ticktype = "detailed", pch = 19, col = NA, border = "gray", nticks = 4)
for (ix in seq(50, 100, 10)) lines (trans3d(x = ix, y = c(-2.5, 12.5), z= 40, pmat = pmat), col = "black") 
for (iy in seq(0, 10, 5)) lines (trans3d(x = c(46, 100), y = iy, z= 40, pmat = pmat), col = "black")

points(trans3d(Coral.cover, Coral.richness, rep(40, 12), pmat = pmat), col = "gray", pch = 19)
points(trans3d(Coral.cover, Coral.richness, fish.rich, pmat = pmat), col = "blue")

xy0 <- trans3d(Coral.cover, Coral.richness, fit$fitted.values, pmat = pmat)
xy1 <- trans3d(Coral.cover, Coral.richness, fish.rich, pmat = pmat)
arrows(xy0[[1]], xy0[[2]], xy1[[1]], xy1[[2]], col = "red", lty = 2, length = 0.1)

Version: 08 Juni, 2020 12:59