data cleaning
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.0 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(ggplot2)
library(psych)##
## Attaching package: 'psych'
##
## The following objects are masked from 'package:ggplot2':
##
## %+%, alphalibrary(car)## Loading required package: carData
##
## Attaching package: 'car'
##
## The following object is masked from 'package:psych':
##
## logit
##
## The following object is masked from 'package:dplyr':
##
## recode
##
## The following object is masked from 'package:purrr':
##
## somelibrary(multcomp)## Loading required package: mvtnorm
## Loading required package: survival
## Loading required package: TH.data
## Loading required package: MASS
##
## Attaching package: 'MASS'
##
## The following object is masked from 'package:dplyr':
##
## select
##
##
## Attaching package: 'TH.data'
##
## The following object is masked from 'package:MASS':
##
## geyserlibrary(cowplot)##
## Attaching package: 'cowplot'
##
## The following object is masked from 'package:lubridate':
##
## stamp#data cleaning
dataset <- read_csv("dataset_123.csv")## Rows: 300 Columns: 16
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): Treatment
## dbl (13): Lizard, year, Aggression, Boldness, sex, avgmass, trialnum, ticks...
## date (2): date, date_infected
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.cleandata<- dataset %>%
dplyr::select(Aggression, avgmass, Treatment, avgtotal) %>% #select all columns included in analysis
drop_na(avgmass, Treatment, avgtotal) # remove all rows with NAs
#creating categories for aggression to simplify viewing/interpretation
cut_points <- c(0, 3, 7, 11)
labels <- c("low", "medium", "high")
# Convert the Aggression column to a categorical variable
cleandata$AggressionCategory <- cut(cleandata$Aggression, breaks = cut_points, labels = labels)testing normality
#RAW AND TRANSFORMED DATA FAILED TO BE NORMAL
mass_hist<- ggplot(data= cleandata, aes(x = avgmass)) + #checking normality visually of avg mass of lizard
geom_histogram(fill = "darkgreen", border = "grey")+ #picking colors- green bc lizards!
theme_bw()+ #simple theme
labs(x= "Average Mass", y= "Frequency")+ #axis lables
theme((axis.title= element_text(face="bold", size = 12))) #font size## Warning in geom_histogram(fill = "darkgreen", border = "grey"): Ignoring
## unknown parameters: `border`mass_hist #looks fairly normal but has some dips## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# qqPlot of y-variable
qqPlot(cleandata$avgmass) #almost all points fall in boundary line!
## [1] 255 295# shapiro-wilks test of y-variable
shapiro.test(cleandata$avgmass) #P is less than 0.05 (alpha), must reject the null that the data is normally distributed##
## Shapiro-Wilk normality test
##
## data: cleandata$avgmass
## W = 0.98855, p-value = 0.01819######testing transformations to attempt to get normal distribution
min(cleandata$avgmass) # not zeros! good to take log## [1] 445cleandata$log_avgmass <- log(cleandata$avgmass) #log transformation
## Transformation 1 - LOG code
logmass_hist<- ggplot(data= cleandata, aes(x = log_avgmass)) +
geom_histogram(fill = "darkgreen", border = "grey")+ #picking colors- green bc lizards!
theme_bw()+ #simple theme
labs(x= "Log Average Mass", y= "Frequency")+ #axis lables
theme((axis.title= element_text(face="bold", size = 12))) #font size## Warning in geom_histogram(fill = "darkgreen", border = "grey"): Ignoring
## unknown parameters: `border`logmass_hist #looks like data is now skewed right... less normal than before## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
qqPlot(cleandata$log_avgmass)#definitely some points out of bounds
## [1] 255 295shapiro.test(cleandata$log_avgmass) #NOPE! still not normal p<0.05##
## Shapiro-Wilk normality test
##
## data: cleandata$log_avgmass
## W = 0.96993, p-value = 6.507e-06## Transformation 1 code
cleandata$sqrt_avgmass <- sqrt(cleandata$avgmass) #sqrt transformation
sqrtmass_hist<- ggplot(data= cleandata, aes(x = sqrt_avgmass)) + #checking normality visually of avg mass of lizard
geom_histogram(fill = "darkgreen", border = "grey")+ #picking colors- green bc lizards!
theme_bw()+ #simple theme
labs(x= "Square Root Average Mass", y= "Frequency")+ #axis lables
theme((axis.title= element_text(face="bold", size = 12))) #font size## Warning in geom_histogram(fill = "darkgreen", border = "grey"): Ignoring
## unknown parameters: `border`sqrtmass_hist #also looks like data is now skewed right... less normal than before## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
qqPlot(cleandata$sqrt_avgmass)#definitely some points out of bounds
## [1] 255 295shapiro.test(cleandata$sqrt_avgmass) #NOPE! still not normal p<0.05##
## Shapiro-Wilk normality test
##
## data: cleandata$sqrt_avgmass
## W = 0.98102, p-value = 0.0005206plots of variable of interest
## plot y~x(numeric)
numscatter <-ggplot(data= cleandata, aes(x= avgtotal, y = avgmass)) +
geom_point()+
theme_bw()+
labs(x= "Average Total Ticks", y= "Average Mass")+
theme(axis.title= element_text(face="bold", size = 12))
numscatter
###plot y~x(categorical)
agro_boxplot<-ggplot(data= cleandata, aes(x= AggressionCategory, y = avgmass)) +
geom_boxplot(fill = "#7bbfc9")+
aes(group = Aggression)+
theme_bw()+
labs(x= "Aggression Score", y= "Average Mass")+
theme(axis.title= element_text(face="bold", size = 12))
agro_boxplot
###plot y~x(categorical)
treat_scatter<-ggplot(data= cleandata, aes(x= Treatment, y = avgmass)) +
geom_boxplot(fill = "#7bbfc9")+
theme_bw()+
labs(x= "Treatment", y= "Average Mass")+
theme(axis.title= element_text(face="bold", size = 12))
treat_scatter
stat test 1- Mann-Whitney-U
#test research question 1- Do lizards have significantly different body masses when infested vs. not?
#see earlier that data does not pass assumptions for parametric test, comparing two different groups of untransformed data
##COMPARING HUMAN INTERFERENCE INFESTION TREATMENT
selected_treatment <- cleandata%>% #only comparing those that were exposed to same conditions
filter(Treatment %in% c("infestation", "no_infestation"))
wilcox.test(selected_treatment$avgmass ~ selected_treatment$Treatment) #nonparametric eq of two-sample t-test, test if medians dif from eachother ##
## Wilcoxon rank sum test with continuity correction
##
## data: selected_treatment$avgmass by selected_treatment$Treatment
## W = 2316, p-value = 0.8112
## alternative hypothesis: true location shift is not equal to 0#COMPARING WILD INFESTATION
selected_treatment2 <- cleandata%>% #only comparing those that were exposed to same conditions
filter(Treatment %in% c("no_treatment"))
cut_points <- c(-Inf, 1, +Inf)
labels <- c("no_infestation", "infestation")
# Convert the Aggression column to a categorical variable
selected_treatment2$notreatmentinfestationCategory <- cut(selected_treatment2$avgtotal, breaks = cut_points, labels = labels)
wilcox.test(selected_treatment2$avgmass ~ selected_treatment2$notreatmentinfestationCategory) #nonparametric eq of two-sample t-test, test if medians dif from eachother ##
## Wilcoxon rank sum test with continuity correction
##
## data: selected_treatment2$avgmass by selected_treatment2$notreatmentinfestationCategory
## W = 2176.5, p-value = 0.05243
## alternative hypothesis: true location shift is not equal to 0#p-value = 0.05243 more than 0.05 fail to rejectstat test 2- correlation
library(psych)
#selecting numerical columns to perform test on
cleandata_numeric<- cleandata%>%
dplyr::select(avgmass, avgtotal) #variables of interest
#check assumptions
pairs.panels(cleandata_numeric, #avg total has hard left sqew
density = TRUE,
cor= TRUE,
lm= TRUE)
#transformation trial
#cleandata$log_avgtotal <- log(cleandata$avgtotal+1) #log transformation
#cleandata$sqrt_avgtotal <- sqrt(cleandata$avgtotal) #log transformation
#take 2 pairs panel
#cleandata_numeric2<- cleandata%>%
# dplyr::select(log_avgmass, log_avgtotal, sqrt_avgmass, sqrt_avgtotal) #variables of interest
#pairs.panels(cleandata_numeric2, #log transform looks best- sticking with that ... though the correlation is very low
# density = TRUE,
# cor= TRUE,
# lm= TRUE)
#str(cleandata_numeric2)
#Run a Pearson's correlation test bc match assumptions with log transformation
cor.test(cleandata$avgmass, #one numeric
cleandata$avgtotal, # second num
method = "spearman", #MUST say pearson or spearman- sp bc not normal
alternative = "two.sided") # bc value can be above OR below threshold ---> 2 sided## Warning in cor.test.default(cleandata$avgmass, cleandata$avgtotal, method =
## "spearman", : Cannot compute exact p-value with ties##
## Spearman's rank correlation rho
##
## data: cleandata$avgmass and cleandata$avgtotal
## S = 3866510, p-value = 0.01468
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.1407659#p-value = 0.01468, reject null!linear regression - do I need to include thise since pairs panel showed no corr?
#linear regression model??
model<- lm(avgmass ~ avgtotal,
data = cleandata)
par(mfrow = c(2,2))
plot(model)
test <- sample(1:nrow(cleandata), 15, replace = F) # pick 15 random rows from air_quality, don't replace
mass_train <- cleandata[-test,] # leave those rows out of training data
mass_test <- cleandata[test,] # use those rows to create a set of test data
## Step 10. Use the best fit model on TRAINING set
fit_mass <- lm(avgmass ~ avgtotal,
data = mass_train)
## Step 11. Use the fitted model with training data to make predictions of your Y outcome
fit_mass_train <- predict(fit_mass,
mass_test)
## Step 12. Visualize how your best fit model did with predictions
plot(mass_test$avgmass, pch = 1, ylab = "Avg Mass")+ # plot actual test data values
points(fit_mass_train, pch = 20, col = "red") # plot the model predictions for those points## integer(0)#clearly, log total avg total is a poor predictor of avg mass