library(tidyverse)
## ── Attaching packages ──────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.0     ✓ purrr   0.3.3
## ✓ tibble  2.1.3     ✓ dplyr   0.8.4
## ✓ tidyr   1.0.2     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.5.0
## ── Conflicts ─────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(ggplot2)
library(plotly)
## Warning: package 'plotly' was built under R version 3.6.2
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout

Get bird fragmentation data

bird_frag <- read.csv(file = "https://mlammens.github.io/ENS-623-Research-Stats/data/Logan_Examples/Chapter9/Data/loyn.csv")
summary(bird_frag)
##      ABUND            AREA            YR.ISOL          DIST       
##  Min.   : 1.50   Min.   :   0.10   Min.   :1890   Min.   :  26.0  
##  1st Qu.:12.40   1st Qu.:   2.00   1st Qu.:1928   1st Qu.:  93.0  
##  Median :21.05   Median :   7.50   Median :1962   Median : 234.0  
##  Mean   :19.51   Mean   :  69.27   Mean   :1950   Mean   : 240.4  
##  3rd Qu.:28.30   3rd Qu.:  29.75   3rd Qu.:1966   3rd Qu.: 333.2  
##  Max.   :39.60   Max.   :1771.00   Max.   :1976   Max.   :1427.0  
##      LDIST            GRAZE            ALT       
##  Min.   :  26.0   Min.   :1.000   Min.   : 60.0  
##  1st Qu.: 158.2   1st Qu.:2.000   1st Qu.:120.0  
##  Median : 338.5   Median :3.000   Median :140.0  
##  Mean   : 733.3   Mean   :2.982   Mean   :146.2  
##  3rd Qu.: 913.8   3rd Qu.:4.000   3rd Qu.:182.5  
##  Max.   :4426.0   Max.   :5.000   Max.   :260.0

Log transform AREA

bird_frag$LOG_AREA <- log10(bird_frag$AREA)

Perform a linear regression

bird_frag_lm <- lm(data = bird_frag, ABUND ~ log10(AREA) + YR.ISOL + log10(DIST) + log10(LDIST) + GRAZE + ALT)

Get model summary

summary(bird_frag_lm)
## 
## Call:
## lm(formula = ABUND ~ log10(AREA) + YR.ISOL + log10(DIST) + log10(LDIST) + 
##     GRAZE + ALT, data = bird_frag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.6506  -2.9390   0.5289   2.5353  15.2842 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -125.69725   91.69228  -1.371   0.1767    
## log10(AREA)     7.47023    1.46489   5.099 5.49e-06 ***
## YR.ISOL         0.07387    0.04520   1.634   0.1086    
## log10(DIST)    -0.90696    2.67572  -0.339   0.7361    
## log10(LDIST)   -0.64842    2.12270  -0.305   0.7613    
## GRAZE          -1.66774    0.92993  -1.793   0.0791 .  
## ALT             0.01951    0.02396   0.814   0.4195    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.384 on 49 degrees of freedom
## Multiple R-squared:  0.6849, Adjusted R-squared:  0.6464 
## F-statistic: 17.75 on 6 and 49 DF,  p-value: 8.443e-11

Look at diagnostic plots

plot(bird_frag_lm)

Steps to make plotly figures to visualize 3-dimensions

x_values <- bird_frag$LOG_AREA %>% 
  round(2)
y_values <- bird_frag$GRAZE %>% 
  round(2)
z_values <- bird_frag$ABUND %>% 
  round(2)
# Define regression plane -------------------------------------------------
# Construct x and y grid elements
x_grid <- seq(from = min(bird_frag$LOG_AREA), to = max(bird_frag$LOG_AREA), length = 50)
y_grid <- seq(from = min(bird_frag$GRAZE), to = max(bird_frag$GRAZE), length = 50)
# Construct z grid by computing
# 1) fitted beta coefficients
# 2) fitted values of outer product of x_grid and y_grid
# 3) extracting z_grid (matrix needs to be of specific dimensions)
beta_hat <- bird_frag %>% 
  #lm(ABUND ~ LOG_AREA + GRAZE, data = .) %>% 
  lm(ABUND ~ LOG_AREA + GRAZE + log10(DIST) + log10(LDIST) + YR.ISOL + ALT, data = .) %>%
  coef()

fitted_values <- crossing(y_grid, x_grid) %>% 
  #mutate(z_grid = beta_hat[1] + beta_hat[2]*x_grid + beta_hat[3]*y_grid)
  mutate(z_grid = beta_hat[1] + beta_hat[2]*x_grid + beta_hat[3]*y_grid + 
           beta_hat[4]*mean(log10(bird_frag$DIST)) + beta_hat[5]*mean(log10(bird_frag$LDIST)) +
           beta_hat[6]*mean(bird_frag$YR.ISOL) + beta_hat[7]*mean(bird_frag$ALT))

z_grid <- fitted_values %>% 
  pull(z_grid) %>%
  matrix(nrow = length(x_grid)) %>%
  t()
# Define text element for each point in plane
text_grid <- fitted_values %>% 
  pull(z_grid) %>%
  round(3) %>% 
  as.character() %>% 
  paste("ABUND: ", ., sep = "") %>% 
  matrix(nrow = length(x_grid)) %>%
  t()
# Plot using plotly -------------------------------------------------------
plot_ly() %>%
  # 3D scatterplot:
  add_markers(
    x = x_values,
    y = y_values,
    z = z_values,
    marker = list(size = 5),
    hoverinfo = 'text',
    text = ~paste(
      "ABUND: ", z_values, "<br>",
      "GRAZE: ", y_values, "<br>",
      "LOG_AREA: ", x_values 
    )
  ) %>%
  # Regression plane:
  add_surface(
    x = x_grid,
    y = y_grid,
    z = z_grid,
    hoverinfo = 'text',
    text = text_grid
  ) %>%
  # Axes labels and title:
  layout(
    title = "3D scatterplot and regression plane",
    scene = list(
      zaxis = list(title = "y: ABUND"),
      yaxis = list(title = "x2: GRAZE"),
      xaxis = list(title = "x1: LOG_AREA")
    )
  )