Explore the difference in rainfall between seeding and non-seeding
experiments by looking at the statistical characteristics (mean, sd,
etc) of the experiments and by using visualization.
clouds <- read.csv("clouds.csv")
library(ggplot2)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
mean_seeded <- mean(clouds$rainfall[clouds$seeding == "yes"])
mean_non_seeded <- mean(clouds$rainfall[clouds$seeding == "no"])
mean_seeded
## [1] 4.634167
mean_non_seeded
## [1] 4.171667
sd_seeded <- sd(clouds$rainfall[clouds$seeding == "yes"])
sd_non_seeded <- sd(clouds$rainfall[clouds$seeding == "no"])
sd_seeded
## [1] 2.776841
sd_non_seeded
## [1] 3.519196
box <- ggplot(clouds, aes(x = seeding, y = rainfall)) +
geom_boxplot() +
labs(x = "Seeding", y = "Rainfall", title = "Distribution of Rainfall by Seeding") +
theme_minimal()
box
Use a t-test to see whether there is a significant difference.
t_test <- t.test(rainfall ~ seeding, data = clouds)
t_test
##
## Welch Two Sample t-test
##
## data: rainfall by seeding
## t = -0.3574, df = 20.871, p-value = 0.7244
## alternative hypothesis: true difference in means between group no and group yes is not equal to 0
## 95 percent confidence interval:
## -3.154691 2.229691
## sample estimates:
## mean in group no mean in group yes
## 4.171667 4.634167
Build a multiple linear regression model to model the effects of
seeding on rainfall.
lm_model <- lm(rainfall ~ seeding + cloudcover + prewetness + echomotion + sne, data = clouds)
summary(lm_model)
##
## Call:
## lm(formula = rainfall ~ seeding + cloudcover + prewetness + echomotion +
## sne, data = clouds)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.1158 -1.7078 -0.2422 1.3368 6.4827
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.37680 2.43432 2.620 0.0174 *
## seedingyes 1.12011 1.20725 0.928 0.3658
## cloudcover 0.01821 0.11508 0.158 0.8761
## prewetness 2.55109 2.70090 0.945 0.3574
## echomotionstationary 2.59855 1.54090 1.686 0.1090
## sne -1.27530 0.68015 -1.875 0.0771 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.855 on 18 degrees of freedom
## Multiple R-squared: 0.3403, Adjusted R-squared: 0.157
## F-statistic: 1.857 on 5 and 18 DF, p-value: 0.1524
anova(lm_model)
## Analysis of Variance Table
##
## Response: rainfall
## Df Sum Sq Mean Sq F value Pr(>F)
## seeding 1 1.283 1.2834 0.1575 0.69613
## cloudcover 1 15.738 15.7377 1.9313 0.18157
## prewetness 1 0.003 0.0027 0.0003 0.98557
## echomotion 1 29.985 29.9853 3.6798 0.07108 .
## sne 1 28.649 28.6491 3.5158 0.07711 .
## Residuals 18 146.677 8.1487
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Which variables appear to influence rainfall the most? SNE
As the suitability index (sne) appears to have a relationship with
rainfall, build two new models relating this variable to rainfall, one
for the seeding experiments, one for the non-seeding experiments.
Compare the coefficients for the two models, and produce a figure
showing the two models.
seeded <- subset(clouds, seeding == "yes")
non_seeded <- subset(clouds, seeding == "no")
seed_model <- lm(rainfall ~ sne, data = seeded)
summary(seeded)
## X seeding time sne
## Min. : 2.00 Length:12 Min. : 1.00 Min. :2.050
## 1st Qu.: 8.75 Class :character 1st Qu.:22.50 1st Qu.:2.850
## Median :13.00 Mode :character Median :33.50 Median :3.275
## Mean :12.83 Mean :36.08 Mean :3.330
## 3rd Qu.:18.50 3rd Qu.:56.00 3rd Qu.:4.032
## Max. :23.00 Max. :82.00 Max. :4.400
## cloudcover prewetness echomotion rainfall
## Min. : 2.800 Min. :0.0730 Length:12 Min. : 1.090
## 1st Qu.: 3.250 1st Qu.:0.1237 Class :character 1st Qu.: 2.683
## Median : 4.650 Median :0.1710 Mode :character Median : 4.530
## Mean : 7.683 Mean :0.3179 Mean : 4.634
## 3rd Qu.: 7.025 3rd Qu.:0.2395 3rd Qu.: 5.468
## Max. :37.900 Max. :1.2670 Max. :11.860
non_seeded_model <- lm(rainfall ~ sne, data = non_seeded)
summary(non_seeded)
## X seeding time sne
## Min. : 1.00 Length:12 Min. : 0.00 Min. :1.300
## 1st Qu.: 6.75 Class :character 1st Qu.:15.75 1st Qu.:2.200
## Median :11.00 Mode :character Median :30.00 Median :3.250
## Mean :12.17 Mean :34.58 Mean :3.008
## 3rd Qu.:17.50 3rd Qu.:53.75 3rd Qu.:3.837
## Max. :24.00 Max. :83.00 Max. :4.650
## cloudcover prewetness echomotion rainfall
## Min. : 2.200 Min. :0.0180 Length:12 Min. : 0.280
## 1st Qu.: 4.500 1st Qu.:0.1875 Class :character 1st Qu.: 1.075
## Median : 6.050 Median :0.2555 Mode :character Median : 4.055
## Mean : 6.808 Mean :0.3362 Mean : 4.172
## 3rd Qu.: 8.375 3rd Qu.:0.4300 3rd Qu.: 5.832
## Max. :13.400 Max. :0.7960 Max. :12.850
coef_seeded <- coef(seed_model)
coef_non_seeded <- coef(non_seeded_model)
coef_data <- data.frame(
Type = c("seeding", "Non-Seeding"),
Intercept = c(coef_seeded, coef_non_seeded),
Slope = c(coef_seeded, coef_non_seeded)
)
coef_data
## Type Intercept Slope
## 1 seeding 12.020237 12.020237
## 2 Non-Seeding -2.218039 -2.218039
## 3 seeding 7.319500 7.319500
## 4 Non-Seeding -1.046371 -1.046371
rainfall_seq <- seq(min(clouds$sne), max(clouds$sne), length.out = 100)
pred_seed <- predict(seed_model, newdata = data.frame(sne = rainfall_seq))
pred_nonseed <- predict(non_seeded_model, newdata = data.frame(sne = rainfall_seq))
plot_data <- data.frame(
sne = rainfall_seq,
rainfall_seed = pred_seed,
rainfall_nonseed = pred_nonseed
)
coef_model <- ggplot(plot_data, aes(x = sne)) +
geom_line(aes(y = rainfall_seed, color = "Seeding"), linewidth = 1.5) +
geom_line(aes(y = rainfall_nonseed, color = "Non-Seeding"), linewidth = 1.5) +
labs(x = "Suitability criterion (SNE)", y = "Rainfall", color = "Type") +
scale_color_manual(values = c("Seeding" = "yellow", "Non-Seeding" = "purple")) +
theme_minimal()
coef_model
What does the difference in slope suggest? As the SNE rises,
rainfall decreases