```{r}
library(tidyverse)
library(ggplot2)
library(plotly)
```


Get plant distribution data

```{r}
c3_plants <- read.csv(file = "https://mlammens.github.io/ENS-623-Research-Stats/data/Logan_Examples/Chapter9/Data/paruelo.csv")
summary(c3_plants)
```

Transform data

```{r}
c3_plants$LAT <- as.numeric(scale(c3_plants$LAT, scale = FALSE))
c3_plants$LONG <- as.numeric(scale(c3_plants$LONG, scale = FALSE))
```

Perform a linear regression

```{r}
c3_plants_lm2 <- lm(data = c3_plants, log10(C3 + 0.1) ~ LONG*LAT)
```

Get model summary

```{r}
summary(c3_plants_lm2)
```


Look at diagnostic plots

```{r}
plot(c3_plants_lm2)
```


Steps to make plotly figures to visualize 3-dimensions

```{r}
x_values <- c3_plants$LAT %>% 
  round(2)
y_values <- c3_plants$LONG %>% 
  round(2)
z_values <- log10(c3_plants$C3 + 0.1) %>% 
  round(2)

```



```{r}
# Define regression plane -------------------------------------------------
# Construct x and y grid elements
x_grid <- seq(from = min(c3_plants$LAT), to = max(c3_plants$LAT), length = 50)
y_grid <- seq(from = min(c3_plants$LONG), to = max(c3_plants$LONG), length = 50)

```


```{r}
# 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 <- c3_plants %>% 
  lm(data = ., log10(C3 + 0.1) ~ LAT*LONG) %>%
  coef()

fitted_values <- crossing(y_grid, x_grid) %>% 
  mutate(z_grid = beta_hat[1] + beta_hat[2]*x_grid + beta_hat[3]*y_grid + beta_hat[4]*x_grid*y_grid)
z_grid <- fitted_values %>% 
  pull(z_grid) %>%
  matrix(nrow = length(x_grid)) %>%
  t()

```




```{r}
# Define text element for each point in plane
text_grid <- fitted_values %>% 
  pull(z_grid) %>%
  round(3) %>% 
  as.character() %>% 
  paste("C3: ", ., sep = "") %>% 
  matrix(nrow = length(x_grid)) %>%
  t()

```

```{r}
# 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(
      "C3: ", z_values, "<br>",
      "LAT: ", y_values, "<br>",
      "LONG: ", 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: C3"),
      yaxis = list(title = "x2: LONG"),
      xaxis = list(title = "x1: LAT")
    )
  )
```

