Villages Principal Component Analysis
Katie
05/01/2021
Set Up
Infection Values Data Frame
Prev.Pos.Dogs.T1 = Prevalence of positive Dog Serology from original data at time point 1
Prev.Pos.Dogs.T3 = Prevalence of positive Dog Serology from original data at time point 3
Inc.MST2 = Incidence of positive Human MST from original data
Mean.Max.Saliva = I calculated the maximum saliva test for each dog and found the average of these maximum scores from each village
Mean.Saliva = I averaged the saliva scores for all tests taken from each village
Prev = The saliva dataset also had serology tests for the dogs so I calculated the prevalence from this
Saliva.Prevlence = If the saliva score was above 916 it counts as positive and below 916 counted as negative. I used all the tests to calculate the saliva prevelence in each village
InfectionValuesdf <- read.csv('InfectionValuesdf.csv')
InfectionValuesdf <- InfectionValuesdf[,1:8]
InfectionValuesdf <- subset(InfectionValuesdf, !is.na(InfectionValuesdf$Mean.Saliva))
row.names(InfectionValuesdf) <- InfectionValuesdf$Group
InfectionValuesdf <- InfectionValuesdf[-1]
InfectionValuesdf## Prev.Pos.Dogs.T1 Prev.Pos.Dogs.T3 Inc.MST2 Mean.Max.Saliva Mean.Saliva
## BCVN 0.7777778 0.9333333 0.06206897 12188.260 7068.542
## BMCA 0.9090909 0.6842105 0.04395604 18520.143 10038.423
## BVSV 0.5833333 0.7250000 0.03900709 6626.440 4106.023
## MGROS 0.5555556 0.8400000 0.01276596 2543.100 2543.100
## MQBA 0.5384615 0.6304348 0.06845238 9086.183 4540.988
## PABE 0.8000000 0.8421053 0.05128205 17818.580 12414.814
## PD 0.7894737 0.7000000 0.05442177 9657.033 6450.597
## PF 0.6666667 0.6666667 0.04687500 9708.000 5116.465
## Prev Saliva.Prevlence
## BCVN 0.7368421 0.7600000
## BMCA 0.6153846 0.7647059
## BVSV 0.4545455 0.6875000
## MGROS 0.0000000 0.1111111
## MQBA 0.5588235 0.8947368
## PABE 0.5714286 0.6666667
## PD 0.5714286 0.6250000
## PF 0.7058824 0.7083333Environmental Variables Data Frame
EnviroVardf <- read.csv('EnviroVardf.csv')
row.names(EnviroVardf) <- EnviroVardf$Group
EnviroVardf <- EnviroVardf[-1]
EnviroVardf## Annual.Precipitation Wettest Driest Seasonal Tmin Tmax Precip
## ABCHAC 2545 465 23 74 233.6667 311.8333 212.0833
## BCVN 2712 513 19 79 236.9167 311.1667 226.0000
## BMCA 2565 472 21 75 234.1667 311.8333 213.7500
## BVSV 2675 502 22 78 236.1667 311.4167 222.9167
## CD 2533 463 21 74 233.8333 312.0000 211.0833
## CL 2804 536 16 80 238.0833 311.0833 233.6667
## CM 2506 457 21 75 234.0833 312.4167 208.8333
## FDR 2484 448 23 73 233.0000 312.4167 207.0000
## JUCACH 2492 451 22 73 233.0833 312.2500 207.6667
## MARU 2549 467 22 75 234.0833 311.9167 212.4167
## MGROS 2548 469 20 75 235.0000 312.0000 212.3333
## MQBA 2751 522 17 79 237.5833 311.0833 229.2500
## PABE 2608 484 21 76 235.5000 311.9167 217.3333
## PD 2595 480 22 76 234.4167 311.5000 216.2500
## PF 2742 520 17 79 237.2500 311.1667 228.5000
## SM 2506 455 21 74 233.3333 311.8333 208.8333
## VC 2634 491 23 77 235.0000 311.3333 219.5000
## VUPD 2586 477 22 75 234.4167 311.6667 215.5000
## NDVIMean NDVIMedian NDVIMax NDVIMin LSTMean LSTMedian LSTMax
## ABCHAC 71171480 72124656 83807144 39128224 300.5819 300.8271 303.0098
## BCVN 68978590 74344300 88563360 23014351 301.3609 301.3182 302.7504
## BMCA 76365818 81642252 89359640 22570000 301.6482 301.5077 303.9783
## BVSV 66396028 71211897 80306307 35999860 301.2537 301.3463 302.7621
## CD 67193928 70285952 85833000 33270751 301.3924 301.0569 303.1886
## CL 61414923 63564123 77102720 21825432 300.6765 300.5820 302.5384
## CM 71262930 75586501 85088944 21250836 302.0870 302.1326 303.5665
## FDR 59531105 67846856 83440456 10942584 299.9074 299.8832 301.1968
## JUCACH 76857983 83674620 90319288 11576548 300.0859 300.0298 302.0617
## MARU 68922772 76347064 87785160 14878260 301.0485 301.0889 302.8358
## MGROS 59298448 64433116 78529904 12650532 303.4854 303.2383 305.1004
## MQBA 68256453 73206280 82936463 31080920 301.2487 301.2910 303.3614
## PABE 57434785 55752292 83119800 24896600 302.8574 302.5098 305.1974
## PD 74303323 76214741 83793944 43079578 300.5899 301.1789 302.7252
## PF 72026454 75880300 86175296 34324358 301.3940 301.3448 303.0893
## SM 67031112 73060000 83226112 13844195 300.1174 300.0407 301.8558
## VC 74922352 76383299 83954072 51465236 301.0097 300.9769 303.6098
## VUPD 68107798 70802201 79956392 17223452 301.1254 301.6835 303.3354
## LSTMin EvapoMean EvapoMedian EvapoMax EvapoMin
## ABCHAC 296.1200 31.48407 30.30290 50.20264 14.176617
## BCVN 299.6700 34.04578 32.61080 59.02321 13.396488
## BMCA 300.3395 40.50005 39.14014 71.94923 15.479762
## BVSV 300.0952 32.44249 29.98976 51.49193 10.128520
## CD 299.6663 28.59371 27.74046 52.54673 8.543652
## CL 298.4493 1479.88998 1479.22331 1488.55842 1472.834867
## CM 300.4409 34.06431 31.96171 54.65224 13.881584
## FDR 299.0244 618.75475 619.07491 632.48888 597.756224
## JUCACH 298.6659 35.28196 36.10414 51.93369 11.716193
## MARU 298.0121 38.78923 38.19121 67.12448 15.347358
## MGROS 301.9213 28.18515 24.99208 49.68626 11.888706
## MQBA 299.2475 255.44454 255.16791 277.84252 236.315641
## PABE 301.4247 26.73978 24.50969 47.55969 11.947536
## PD 293.0273 38.08400 36.30581 58.62634 14.902667
## PF 299.5918 33.13820 31.46552 59.30126 8.860002
## SM 298.2349 3274.23554 3274.23863 3274.24589 3274.221180
## VC 299.6946 33.92231 33.79157 53.97996 13.031496
## VUPD 295.1025 30.39886 31.15986 49.19620 14.207052PCA for Infection Values
pca <- prcomp(InfectionValuesdf , scale = TRUE)
#plot(pca$x[,1], pca$x[,2])
pca.var <- pca$sdev^2
pca.var.per <- round(pca.var/sum(pca.var)*100, 1)Scree Plot for Infection Values PCA
barplot(pca.var.per, main="Scree Plot", xlab="Principal Component", ylab="Percent Variation")
Results for infection values PCA
pca.data <- data.frame(Sample=rownames(pca$x),
X=pca$x[,1],
Y=pca$x[,2])
ggplot(data=pca.data, aes(x=X, y=Y, label=Sample)) +
geom_text() +
xlab(paste("PC1 - ", pca.var.per[1], "%", sep="")) +
ylab(paste("PC2 - ", pca.var.per[2], "%", sep="")) +
theme_bw() +
ggtitle("PCA for Villages")
Top 5 infection values
loading_scores <- pca$rotation[,1]
infection_scores <- abs(loading_scores)
infection_score_ranked <- sort(infection_scores, decreasing=TRUE)
top_5_infection_values <- names(infection_score_ranked[1:5])
top_5_infection_values ## [1] "Mean.Max.Saliva" "Prev" "Saliva.Prevlence" "Mean.Saliva"
## [5] "Inc.MST2"PCA for Environmental Values
pca <- prcomp(EnviroVardf , scale = TRUE)
#plot(pca$x[,1], pca$x[,2])
pca.var <- pca$sdev^2
pca.var.per <- round(pca.var/sum(pca.var)*100, 1)Scree Plot for Environmental Values PCA
barplot(pca.var.per, main="Scree Plot", xlab="Principal Component", ylab="Percent Variation")
Results for infection values PCA
pca.data <- data.frame(Sample=rownames(pca$x),
X=pca$x[,1],
Y=pca$x[,2])
ggplot(data=pca.data, aes(x=X, y=Y, label=Sample)) +
geom_text() +
xlab(paste("PC1 - ", pca.var.per[1], "%", sep="")) +
ylab(paste("PC2 - ", pca.var.per[2], "%", sep="")) +
theme_bw() +
ggtitle("PCA for Villages")
Top 5 environmental values
loading_scores <- pca$rotation[,1]
environmental_scores <- abs(loading_scores)
environmental_score_ranked <- sort(environmental_scores, decreasing=TRUE)
top_5_environmental_values <- names(environmental_score_ranked[1:5])
top_5_environmental_values ## [1] "Tmin" "Wettest" "Annual.Precipitation"
## [4] "Precip" "Seasonal"```