# -*- coding: utf-8 -*-
"""passnyc-numeric-and-categorical-variables-r.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1J2ntDrNzoCf2xp0eIwKDgOFUWMgVQqqZ

<h1 class='font-effect-3d' style='color:#37c9e1; font-family:Orbitron;'> &#x1F310; &nbsp; Code Library, Styling, and Links</h1>
<details><summary style='color:#37c9e1; font-family:Orbitron;'>Github</summary><br/>

The current notebooks
    
&#x1F4D8; &nbsp; [Python Version](https://github.com/OlgaBelitskaya/kaggle_notebooks/blob/master/passnyc-numeric-and-categorical-variables.ipynb)

&#x1F4D8; &nbsp; [R Version](https://github.com/OlgaBelitskaya/kaggle_notebooks/blob/master/passnyc-numeric-and-categorical-variables-r.ipynb)
    
The previous notebooks

&#x1F4D8; &nbsp; [Python Version](https://github.com/OlgaBelitskaya/kaggle_notebooks/blob/master/passnyc-data-exploration.ipynb)

&#x1F4D8; &nbsp; [R Version](https://github.com/OlgaBelitskaya/kaggle_notebooks/blob/master/passnyc-data-exploration-r.ipynb)

<details><summary style='color:#37c9e1; font-family:Orbitron;'>Useful Links</summary><br/>

&#x1F4E1; &nbsp; [School Quality Reports. Educator Guide](http://schools.nyc.gov/NR/rdonlyres/967E0EE1-7E5D-4E47-BC21-573FEEE23AE2/0/201516EducatorGuideHS9252017.pdf)
    
&#x1F4E1; &nbsp; [New York City Department of Education](https://www.schools.nyc.gov)

&#x1F4E1; &nbsp; [NYC OpenData](https://opendata.cityofnewyork.us/)

&#x1F4E1; &nbsp; [Pandas Visualization](https://pandas.pydata.org/pandas-docs/stable/visualization.html)
    
&#x1F4E1; &nbsp; [Pandas Styling](https://pandas.pydata.org/pandas-docs/stable/style.html)

&#x1F4E1; &nbsp; [ggplot2](https://ggplot2.tidyverse.org/)

&#x1F4E1; &nbsp; [R Tutorial](https://www.tutorialspoint.com/r/index.htm)
"""

library(IRdisplay); library(repr)
library(tidyverse); library(zoo)
library(grid); library(gridExtra)
library(sf); library(Rtsne)
library(corrplot)
cols<-c("#3636ff","#ff3636","#6f6fff","#ff6f6f")

display_html("<style> 
@import url('https://fonts.googleapis.com/css?family=Orbitron|Roboto&effect=3d');
body {background-color: gainsboro;} 
h4 {color:#818286; font-family:Roboto;}
span {color:black; text-shadow:4px 4px 4px #aaa;}
div.output_prompt,div.output_area pre {color:slategray;}
div.input_prompt,div.output_subarea {color:#37c9e1;}      
div.output_stderr pre {background-color:gainsboro;}  
div.output_stderr {background-color:slategrey;}
</style>")

b <- "#3636ff"
r <- "#ff3636"
lr <- "#ff6f6f"
lb <- "#6f6fff"

"""<h1 class='font-effect-3d' style='color:#37c9e1; font-family:Orbitron;'> &#x1F310; &nbsp; Data Loading & Preprocessing<h1>"""

dir1<-"../input/data-science-for-good/"
dir2<-"../input/ny-school-districts/"
file1<-"2016 School Explorer.csv"
file2<-"D5 SHSAT Registrations and Testers.csv"
school_explorer<-read.csv(paste(dir1,file1,sep=""))
d5_shsat<-read.csv(paste(dir1,file2,sep=""))
school_explorer<-school_explorer[,4:161]
c(paste0('school_explorer =>',list(dim(school_explorer))),
  paste0('d5_shsat =>',list(dim(d5_shsat))))

cat(colnames(school_explorer)[1:38],sep=", ")

school_explorer[school_explorer=="N/A"]<-NA
school_explorer$School.Name<-
as.character(school_explorer$School.Name)
school_explorer$School.Name[c(428,1024,713,909)]<- 
c('P.S. 212 D12','P.S. 212 D30','P.S. 253 D21','P.S. 253 D27')
school_explorer$School.Income.Estimate<-
as.character(school_explorer$School.Income.Estimate)
school_explorer$School.Income.Estimate<-
sub("\\$","",school_explorer$School.Income.Estimate)
school_explorer$School.Income.Estimate<-
sub(",","",school_explorer$School.Income.Estimate)
school_explorer$School.Income.Estimate<-
as.numeric(school_explorer$School.Income.Estimate)
school_explorer$School.Income.Estimate<-
na.approx(school_explorer$School.Income.Estimate)
percent_list<-c('Percent.ELL',
                'Percent.Asian','Percent.Black','Percent.Hispanic',
                'Percent.Black...Hispanic','Percent.White',
                'Student.Attendance.Rate',
                'Percent.of.Students.Chronically.Absent',
                'Rigorous.Instruction..','Collaborative.Teachers..',
                'Supportive.Environment..','Effective.School.Leadership..',
                'Strong.Family.Community.Ties..','Trust..')
target_list<-c("Average.ELA.Proficiency","Average.Math.Proficiency")
economic_list<-c("Economic.Need.Index","School.Income.Estimate")
rating_list<-c('Rigorous.Instruction.Rating',
               'Collaborative.Teachers.Rating',
               'Supportive.Environment.Rating',
               'Effective.School.Leadership.Rating',
               'Strong.Family.Community.Ties.Rating','Trust.Rating',
               'Student.Achievement.Rating')
nastr2num<-function(x) {x<-as.numeric(sub("%","",x))}
for (el in c(percent_list,target_list,"Economic.Need.Index")) {
    school_explorer[el]<-sapply(school_explorer[el],nastr2num)
    school_explorer[el]<-na.approx(school_explorer[el]) }
for (el in rating_list) {
    vtable<-data.frame(table(school_explorer[el]))
    mvalue<-as.character(vtable$Var1[vtable$Freq==max(vtable$Freq)])
    school_explorer[el][is.na(school_explorer[el])]<-mvalue
    school_explorer[,el]<-factor(school_explorer[,el]) }

sum(is.na(school_explorer))

fl<-c('District',target_list,economic_list)
districts<-
data.frame(aggregate(.~District,school_explorer[fl],mean))
districts

geo_districts<-sf::st_read(paste(dir2,"nysd.shp",sep=""))
geo_districts$centroid<-st_centroid(geo_districts$geometry)
geo_districts<-st_transform(geo_districts,crs=4326)
geo_districts$centroid<-st_transform(geo_districts$centroid,crs=4326)
geo_districts<-data.frame(geo_districts)
geo_districts<-geo_districts[order(geo_districts$SchoolDist),]
rownames(geo_districts)<-c(1:33)
coords<-do.call(rbind,st_geometry(geo_districts$centroid)) %>% 
            as_tibble() %>% 
            setNames(c("lon","lat"))
geo_districts[c('lon','lat')]<-coords
head(geo_districts,11)

options(repr.plot.width=11,repr.plot.height=8)
ggplot(geo_districts,aes(fill=c(1:33)))+ 
    geom_sf(data=geo_districts,size=.5,
            color="slategray",alpha=.8)+ 
    geom_label(aes(x=lon,y=lat,label=SchoolDist),size=3)+
    scale_fill_distiller(name="Districts",palette="Spectral")+
    ggtitle("School Districts")+theme_bw()

"""<h1 class='font-effect-3d' style='color:#37c9e1; font-family:Orbitron;'> &#x1F310; &nbsp; Classification of Variables</h1>"""

cat('Numeric Variables')
numeric_list<-unlist(lapply(school_explorer,is.numeric))
numeric_list<-colnames(school_explorer[numeric_list])
numeric_list[1:25]

summary(school_explorer[numeric_list][,1:9])

cat('Categorical Variables')
category_list<-!unlist(lapply(school_explorer,is.numeric))
category_list<-colnames(school_explorer[category_list])
category_list

summary(school_explorer[category_list][,1:9])

"""Of course, the variables `SED Code`, `District`, `Zip` are categorical. 

Just categories are denoted by numeric values.

Let's convert string values of other categorical features into numeric.
"""

numeric_list<-numeric_list[c(6:21,24:143)]
numeric_list[1:20]

display_html('<p style="color:#3636ff; font-family:Roboto; font-size:150%; text-shadow:5px 5px 5px #aaa;">
Unique Values of Categorical Variables</p>')
cat('District')
unique(school_explorer$District)
cat('City')
unique(school_explorer$City)
cat('Grades')
unique(school_explorer$Grades)
cat('Community School')
unique(school_explorer$Community.School.)
cat('Trust Rating')
unique(school_explorer$Trust.Rating)

category_list<-c('District','City','Grades',
                 'Community.School.',rating_list)
fac2num<-function(x) {as.numeric(as.factor(x))}
school_explorer[,category_list]<-
lapply(school_explorer[,category_list],fac2num)
display_html('<p style="color:#3636ff; font-family:Roboto; font-size:150%; text-shadow:5px 5px 5px #aaa;">
Unique Values of Categorical Variables after Transformation</p>')
cat('District')
unique(school_explorer$District)
cat('City')
unique(school_explorer$City)
cat('Grades')
unique(school_explorer$Grades)
cat('Community School')
unique(school_explorer$Community.School.)
cat('Trust Rating')
unique(school_explorer$Trust.Rating)

"""<h1 class='font-effect-3d' style='color:#37c9e1; font-family:Orbitron;'> &#x1F310; &nbsp; Clustering</h1>"""

n<-10; fl<-c(target_list,economic_list)
clusters<-kmeans(school_explorer[,fl],n,iter.max=500,
                 algorithm="Hartigan-Wong")
school_explorer['Education.Clusters']<-as.factor(clusters$cluster)
# by education results and economic indicators
str(clusters)

"""Let's arrange these clusters in accordance with educational achievements."""

fl<-c('Education.Clusters',target_list,economic_list)
ed_clusters<-data.frame(
    aggregate(.~Education.Clusters,school_explorer[fl],mean))
ed_clusters<-ed_clusters[order(ed_clusters$Average.Math.Proficiency),]
rownames(ed_clusters)<-c(1:10)
ed_clusters

cl<-c("01","02","03","04","05",
      "06","07","08","09","10")
mapping<-setNames(cl,c(ed_clusters$Education.Clusters))                    
school_explorer$Education.Clusters<-
mapping[unlist(school_explorer$Education.Clusters)]
ed_clusters$Education.Clusters<-cl
ed_clusters

"""The location of the results with good quality by districts and cluster types of schools can be assessed using visualization."""

options(repr.plot.width=11,repr.plot.height=5)
p1<-ggplot(school_explorer,aes(x=factor(District),
                               fill=Education.Clusters))+
          geom_bar(position="dodge",width=.6)+theme_bw()+
          xlab("District")+ggtitle('Education Clusters by Districts')+
          scale_fill_brewer(name="Education Clusters",palette="Spectral")+
          theme(plot.title=element_text(size=15,hjust=.5),
                axis.text.x=element_text(size=10,angle=0))
p1

options(repr.plot.width=11,repr.plot.height=7)
p2<-ggplot(school_explorer)+
        geom_point(aes(x=Longitude,y=Latitude,color=Education.Clusters))+ 
        scale_color_brewer(name="Education Clusters",palette="Spectral")+
        geom_sf(data=geo_districts,size=.5,color="slategray",alpha=.2)+
        theme_bw()+ggtitle('Education Clusters by Schools')+
        theme(plot.title=element_text(size=15,hjust=.5),
              axis.text.x=element_text(size=10,angle=0))
p2

options(repr.plot.width=5,repr.plot.height=4)
plot_variable1<-function(v) {
    ggplot(school_explorer,
           aes(x=factor(Community.School.),
               fill=factor(Community.School.)))+
    geom_boxplot(aes(y=c(factor(school_explorer[,v]))))+
    ggtitle(gsub("\\."," ",v))+
    theme_bw()+xlab("")+ylab("")+
    scale_fill_manual(name="Community School?",
                      values=c(cols[2],cols[1]),
                      labels=c("No","Yes"))+
    scale_x_discrete(labels=c("No","Yes"))+
    theme(plot.title=element_text(size=15,hjust=.5),
          axis.text.x=element_text(size=10,angle=0)) }
plot_variable1('Education.Clusters')

n<-10; fl<-c(target_list,percent_list)
clusters<-kmeans(school_explorer[,fl],n,iter.max=500,
                 algorithm="Hartigan-Wong")
school_explorer['Education.Clusters']<-as.factor(clusters$cluster)
# by education results, the social and demographic situation
str(clusters)

fl<-c('Education.Clusters',target_list,economic_list)
ed_clusters<-data.frame(
    aggregate(.~Education.Clusters,
              school_explorer[fl],mean))
ed_clusters<-ed_clusters[order(ed_clusters$Average.Math.Proficiency),]
rownames(ed_clusters)<-c(1:10)
ed_clusters

cl<-c("01","02","03","04","05","06","07","08","09","10")
mapping<-setNames(cl,c(ed_clusters$Education.Clusters))            
school_explorer$Education.Clusters<-
mapping[unlist(school_explorer$Education.Clusters)]
ed_clusters$Education.Clusters<-cl
ed_clusters

options(repr.plot.width=11,repr.plot.height=5); p1

options(repr.plot.width=11,repr.plot.height=7); p2

options(repr.plot.width=5,repr.plot.height=4)
plot_variable1('Education.Clusters')

"""The distribution by clusters has certain differences for two cases, but the problem regions are clearly the same.
<h1 class='font-effect-3d' style='color:#37c9e1; font-family:Orbitron;'>&#x1F310; &nbsp; Manifold Learning</h1>
Let's divide schools into educational clusters using numeric and categorical variables.
"""

n<-10; fl<-c(target_list,numeric_list,category_list)
clusters<-kmeans(school_explorer[,fl],n,iter.max=500,
                 algorithm="Hartigan-Wong")
school_explorer['Education.Clusters']<-as.factor(clusters$cluster)
# by preprocessed features
str(clusters)

fl<-c('Education.Clusters',target_list,economic_list)
ed_clusters<-data.frame(
    aggregate(.~Education.Clusters,school_explorer[fl],mean))
ed_clusters<-ed_clusters[order(ed_clusters$Average.Math.Proficiency),]
rownames(ed_clusters)<-c(1:10)
ed_clusters

cl<-c("01","02","03","04","05","06","07","08","09","10")
mapping<-setNames(cl,c(ed_clusters$Education.Clusters))                    
school_explorer$Education.Clusters<-
mapping[unlist(school_explorer$Education.Clusters)]
ed_clusters$Education.Clusters<-cl
ed_clusters

options(repr.plot.width=11,repr.plot.height=5); p1

options(repr.plot.width=11,repr.plot.height=7); p2

"""Now we can check whether the multidimensional data form a certain structure or don't. 

For this purpose, a probability algorithm and a transformation to two-dimensional space were used. 

The color on the chart still marks educational clusters.
"""

fl<-c(target_list,numeric_list,category_list)
tsne<-Rtsne(school_explorer[,fl],dims=2,perplexity=50,
            verbose=TRUE,max_iter=500)
data_tsne<-data.frame(tsne$Y)
data_tsne$cluster<-school_explorer$Education.Clusters

options(repr.plot.width=7,repr.plot.height=5)
ggplot(data=data_tsne,aes(x=X1,y=X2,color=cluster))+ 
    theme_bw()+geom_point()+
    scale_color_brewer(name="Education Clusters",
                       palette="Spectral")+
    ggtitle("TSNE Patterns in the Data")+
    theme(plot.title=element_text(size=15,hjust=.5))

"""It seems like this clustering is done correctly.
<h1 class='font-effect-3d' style='color:#37c9e1; font-family:Orbitron;'> &#x1F310; &nbsp; Feature Correlation</h1>
Many numeric variables demonstrate a strong correlation.
"""

display_html('<p style="color:#3636ff; font-family:Roboto; font-size:150%; text-shadow:5px 5px 5px #aaa;">
Correlation Matrix</p>')
fl<-c(target_list,numeric_list[1:16])
corr_matrix<-cor(school_explorer[,fl],
                 use="complete.obs",method="pearson")
options(repr.plot.width=7,repr.plot.height=7)
corrplot(corr_matrix,method="pie",order="hclust",
         tl.col="black",tl.cex=.7)

"""Education results (ELA & Math) are correlated really strongly so they can be combined into one indicator."""

fl<-c(target_list,numeric_list,category_list)
corr_matrix<-cor(school_explorer[,fl],
                 use="complete.obs",method="pearson")
corr_ela<-corr_matrix[,1]
corr_ela<-data.frame(corr_ela[order(-abs(corr_ela))][2:21])
corr_math<-corr_matrix[,2]
corr_math<-data.frame(corr_math[order(-abs(corr_math))][2:21])
corr_ela; corr_math

"""<h1 class='font-effect-3d' style='color:#37c9e1; font-family:Orbitron;'> &#x1F310; &nbsp; Let's Go Ahead</h1>
To be continued...

&#x1F4D8; &nbsp;  [PASSNYC. Comparing All Districts with 5th District](https://www.kaggle.com/olgabelitskaya/passnyc-comparing-all-districts-with-5th-district/)
"""