anal.r
#setwd('~/uhk/dev/ar_ondrej_plachy/institution/')
data = read.table('data.tsv',sep="\t",header=TRUE)
data$time = data$year - 1983
data$random1 = rnorm(336, mean=0, sd=0.01)
data$random2 = rnorm(336, mean=0, sd=0.01)
attach(data)
fbalm2 =
ifelse(balance_of_parties_mandates<= 0.25, '0.0-0.25',
ifelse(balance_of_parties_mandates <= 0.5, '0.25-0.5',
ifelse(balance_of_parties_mandates <= 0.75, '0.5-0.75',
'0.75-1.0')))
model0a = glm(event~time + provincia + institution, family="binomial")
model1a2 = glm(event~time + provincia + fbalm2 + institution, family="binomial")
summary(model1a2)
#likelihood ratio test
library(epicalc)
lrtest (model0a, model1a2)
model1a = glm(event~time + provincia + balance_of_parties + institution, family="binomial")
model1b = glm(event~time + provincia + effective_number_of_parties + institution, family="binomial")
model1c = glm(event~time + provincia + raes_index + institution, family="binomial")
model2a = glm(event~time + provincia + balance_of_parties + effective_number_of_parties + institution, family="binomial")
model2b = glm(event~time + provincia + balance_of_parties + raes_index + institution, family="binomial")
model3a = glm(event~time + provincia + balance_of_parties + effective_number_of_parties + institution + balance_of_parties*effective_number_of_parties, family="binomial")
model3b = glm(event~time + provincia + balance_of_parties + raes_index + institution + balance_of_parties*effective_number_of_parties, family="binomial")
library("splines")
model4a = glm(event~time + provincia + bs(balance_of_parties) + bs(effective_number_of_parties) + institution + bs(balance_of_parties)*bs(effective_number_of_parties), family="binomial")
#factors by quantiles
#install.packages("gtools")
library("gtools")
fbalm = quantcut(balance_of_parties_mandates, q=c(0,0.25,0.5,0.75,1))
fbalm2 =
ifelse(balance_of_parties_mandates<= 0.2, '0.0-0.2',
ifelse(balance_of_parties_mandates <= 0.4, '0.2-0.4',
ifelse(balance_of_parties_mandates <= 0.6, '0.4-0.6',
ifelse(balance_of_parties_mandates <= 0.8, '0.6-0.8',
'0.8-1.0'))))
fbalm2 =
ifelse(balance_of_parties_mandates<= 0.25, '0.0-0.25',
ifelse(balance_of_parties_mandates <= 0.5, '0.25-0.5',
ifelse(balance_of_parties_mandates <= 0.75, '0.5-0.75',
'0.75-1.0')))
model1a = glm(event~time + provincia + fbalm + institution, family="binomial")
model1a2 = glm(event~time + provincia + fbalm2 + institution, family="binomial")
model0a = glm(event~time + provincia + institution, family="binomial")
# http://www.cookbook-r.com/Graphs/Scatterplots_%28ggplot2%29/
library("ggplot2")
ggplot(data,aes(x=balance_of_parties+random1,y=effective_number_of_parties+random2,color=as.factor(event))) + geom_point(shape=1)
ggplot(data,aes(x=balance_of_parties+random1,y=effective_number_of_parties+random2, color=as.factor(event), size=as.factor(event))) + geom_point(shape=1)
# http://stackoverflow.com/questions/11875941/3d-equivalent-of-the-curve-function-in-r
library(emdbook)
curve3d(-2.4966 -11.69859*x - 2.54043*y + 5.44407*x*y,from=c(0,0), to=c(1,5),sys3d=c("contour"))
curve3d(exp(-2.4966 -11.69859*x - 2.54043*y + 5.44407*x*y),from=c(0,0), to=c(1,5),sys3d=c("contour"))
#library("splines", lib.loc="/usr/lib/R/library")
termplot(model3b)
dev.off()
# contour lines:
zmodel3b = function(x,y) { #balance_of_parties, effective_number_of_parties
return (-2.4966 -11.69859*x - 2.54043*y + 5.44407*x*y)
}
x = seq(0.15, 1.025, by = 0.05)
y = seq(0.5, 4.5, by = 0.2)
z = array(0,c(length(x),length(y)))
i = 1
j = 1
for (iks in x) {
for (yps in y) {
z[i,j] = zmodel3b(iks,yps)
j = j + 1
}
i = i + 1
j = 1
}
a = contourLines(x,y,z)
xy = paste(unlist(lapply(a, function(z) {
xs = paste(round(z$x, 3), collapse = ",")
ys = paste(round(z$y, 3), collapse = ",")
sprintf("{\n \"x\": [%s],\n \"y\": [%s]\n}", xs, ys)
})), collapse = ", \n")
cat("<script>", sprintf("var data = [%s]", xy), "</script>", sep = "\n")
# minus controled
data$event_minus_controled = exp(resid(model3b,type="resp") + (-11.69859*balance_of_parties - 2.54043*effective_number_of_parties + 5.44407*balance_of_parties*effective_number_of_parties))
data$fitted = model3b$fitted
data$fitted_minus_controled = exp(-11.69859*balance_of_parties - 2.54043*effective_number_of_parties + 5.44407*balance_of_parties*effective_number_of_parties)
write.csv(data,file="data.csv")
#likelihood ratio test
library(epicalc)
lrtest (model0a, model1a)
anal_rCharts.r
# this assumes that data is collected into
# data and contourlines
# see script anal.R for the code for data and calculations
# there appears to be two data files
# generated to produce the scatterplot with contour lines
# data.csv contains the data for the points on the scatter plot
# contour_lines.json has the data to render the contour lines
# in the original anal.R contour lines
# came from
#a = contourLines(x,y,z)
#xy = paste(unlist(lapply(a, function(z) {
# xs = paste(round(z$x, 3), collapse = ",")
# ys = paste(round(z$y, 3), collapse = ",")
# sprintf("{\n \"x\": [%s],\n \"y\": [%s]\n}", xs, ys)
#})), collapse = ", \n")
#cat("<script>", sprintf("var data = [%s]", xy), "</script>", sep = "\n")
# here we will recreate the chart with rCharts
# rCharts provides several features that might be useful
# for embedding, sharing, and publishing
# also this allows more robust reusability
# we will use an approach much like Ben Hunter
# in his post http://mostlyconjecture.com/2014/05/03/chord-diagrams-with-rcharts/
# to start let's load the data as if we have completed our analysis and ready to plot
require(jsonlite)
dataContours <- fromJSON("contour_lines.json")
dataPoints <- read.csv(
"data.csv"
,stringsAsFactors = F
)
require(rCharts)
# we will make a ref class subset of rCharts called contourLines
contourLines = setRefClass('contourLines', contains = 'rCharts', methods = list(
initialize = function(){
callSuper()
LIB <<- get_lib(".")
lib <<- "contour_lines"
templates$script <<- '
<script>
function draw{{chartId}}(){
var params = {{{ chartParams }}}
var margin = (params.margin) ? params.margin : {top: 20, right: 20, bottom: 30, left: 50},
width = params.width - margin.left - margin.right,
height = params.height - margin.top - margin.bottom;
var x = d3.scale.linear()
.range([0, width])
.domain([0.15,1.025]);
var y = d3.scale.linear()
.range([height, 0])
.domain([0.5,4.5]);
var xAxis = d3.svg.axis()
.scale(x)
.orient("bottom");
var yAxis = d3.svg.axis()
.scale(y)
.orient("left");
var svg = d3.select("#"+params.id).append("svg")
.attr("width", width + margin.left + margin.right)
.attr("height", height + margin.top + margin.bottom)
.append("g")
.attr("transform", "translate(" + margin.left + "," + margin.top + ")");
svg.append("g")
.attr("class", "x axis")
.attr("transform", "translate(0," + height + ")")
.call(xAxis)
.append("text")
.attr("x", x(1))
.attr("dx", ".71em")
.attr("dy", "-.5em")
.style("text-anchor", "end")
.text("Balance of parties");
svg.append("g")
.attr("class", "y axis")
.call(yAxis)
.append("text")
.attr("transform", "rotate(-90)")
.attr("y", 6)
.attr("dy", ".71em")
.style("text-anchor", "end")
.text("Effective number of parties");
//contour lines based on http://vis.supstat.com/2012/11/contour-plots-with-d3-and-r/
var line = d3.svg.line()
.x(function(d) { return x(d.x); })
.y(function(d) { return y(d.y); })
.interpolate("basis") //https://github.com/mbostock/d3/wiki/SVG-Shapes#line_interpolate
.tension(1);
svg.selectAll(".contour-line")
.data(params.data.cl.map(function(d) {
return d3.range(d.x.length).map(function(i) {
return {x: d.x[i], y: d.y[i]};
});
}))
.enter().append("svg:path")
.attr("d", line)
.attr("class","contour-line")
svg.selectAll(".event")
.data(params.data.points)
.enter().append("circle")
.attr("cx", function(d) {
return x(parseFloat(d[params.x]) + parseFloat(d[params.xjitter])/2)
})
.attr("cy", function(d) {
return y(parseFloat(d[params.y]) + parseFloat(d[params.yjitter])*2)
})
.attr("r", function(d) {
return Math.round(parseFloat(d[params.radius])*100 + 1);
})
.attr("class","point noevent")
.filter(function(d) { return (d[params.event] == 1) })
.attr("class","point event")
.attr("r", function(d) {
return Math.min(30,Math.round(parseFloat(d[params.radius])*100 + 1));
})
}
draw{{chartId}}();
</script>
'
},
getPayload = function(chartId){
list(
chartParams = jsonlite:::toJSON(params, auto_unbox=T)
, chartId = chartId
, lib = basename(lib)
, liburl = LIB$url
)
}
))
clPlot <- contourLines$new()
clPlot$set(
data = list(
cl = dataContours,
points = dataPoints
),
x = "balance_of_parties",
xjitter = "random1",
y = "effective_number_of_parties",
yjitter = "random2",
event = "event",
radius = "event_minus_controled",
height = 300,
width = 500
)
clPlot
#clPlot$save("index.html",cdn=F)
contour_lines.json
[{
"x": [0.15,0.151],
"y": [4.495,4.5]
},
{
"x": [0.917,0.95,0.979,1],
"y": [0.5,0.612,0.7,0.756]
},
{
"x": [0.15,0.18,0.2,0.207,0.229],
"y": [3.915,4.1,4.246,4.3,4.5]
},
{
"x": [0.806,0.85,0.853,0.9,0.914,0.95,1,1],
"y": [0.5,0.69,0.7,0.858,0.9,0.992,1.1,1.1]
},
{
"x": [0.15,0.189,0.2,0.225,0.25,0.252,0.274,0.292,0.3,0.307],
"y": [3.335,3.5,3.557,3.7,3.882,3.9,4.1,4.3,4.402,4.5]
},
{
"x": [0.694,0.7,0.726,0.75,0.767,0.8,0.825,0.85,0.9,0.909,0.95,1],
"y": [0.5,0.54,0.7,0.824,0.9,1.022,1.1,1.169,1.282,1.3,1.372,1.445]
},
{
"x": [0.15,0.2,0.211,0.25,0.265,0.3,0.3,0.325,0.343,0.35,0.357,0.368,0.377,0.385],
"y": [2.755,2.868,2.9,3.034,3.1,3.3,3.3,3.5,3.7,3.794,3.9,4.1,4.3,4.5]
},
{
"x": [0.583,0.599,0.6,0.62,0.65,0.65,0.693,0.7,0.75,0.762,0.8,0.85,0.894,0.9,0.95,1],
"y": [0.5,0.7,0.71,0.9,1.1,1.103,1.3,1.327,1.472,1.5,1.573,1.648,1.7,1.706,1.752,1.789]
},
{
"x": [0.15,0.2,0.25,0.3,0.35,0.4,0.413,0.443,0.45,0.452,0.456,0.458,0.46,0.461,0.461,0.462,0.462,0.463,0.463],
"y": [2.175,2.179,2.186,2.198,2.219,2.271,2.3,2.5,2.638,2.7,2.9,3.1,3.3,3.5,3.7,3.9,4.1,4.3,4.5]
},
{
"x": [0.472,0.472,0.473,0.474,0.476,0.479,0.485,0.499,0.5,0.55,0.6,0.633,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1],
"y": [0.5,0.7,0.9,1.1,1.3,1.5,1.7,1.9,1.905,2.051,2.088,2.1,2.104,2.114,2.12,2.124,2.128,2.13,2.132,2.134]
},
{
"x": [0.36,0.35,0.345,0.326,0.3,0.299,0.26,0.25,0.2,0.196,0.15],
"y": [0.5,0.644,0.7,0.9,1.095,1.1,1.3,1.339,1.491,1.5,1.594]
},
{
"x": [1,0.967,0.95,0.9,0.85,0.8,0.785,0.75,0.7,0.7,0.651,0.65,0.619,0.6,0.597,0.58,0.567,0.557,0.55,0.548,0.541],
"y": [2.478,2.5,2.512,2.554,2.607,2.675,2.7,2.768,2.9,2.901,3.1,3.106,3.3,3.465,3.5,3.7,3.9,4.1,4.255,4.3,4.5]
},
{
"x": [0.249,0.219,0.2,0.179,0.15],
"y": [0.5,0.7,0.802,0.9,1.014]
},
{
"x": [1,0.95,0.945,0.9,0.85,0.844,0.8,0.779,0.75,0.733,0.7,0.698,0.672,0.651,0.65,0.634,0.619],
"y": [2.822,2.892,2.9,2.978,3.086,3.1,3.226,3.3,3.417,3.5,3.688,3.7,3.9,4.1,4.108,4.3,4.5]
},
{
"x": [1,0.95,0.938,0.9,0.868,0.85,0.817,0.8,0.777,0.75,0.745,0.719,0.7,0.698],
"y": [3.167,3.272,3.3,3.402,3.5,3.565,3.7,3.777,3.9,4.065,4.1,4.3,4.475,4.5]
},
{
"x": [1,0.95,0.935,0.9,0.882,0.85,0.839,0.804,0.8,0.776],
"y": [3.511,3.652,3.7,3.826,3.9,4.044,4.1,4.3,4.328,4.5]
},
{
"x": [1,0.986,0.95,0.933,0.9,0.89,0.854],
"y": [3.856,3.9,4.032,4.1,4.249,4.3,4.5]
},
{
"x": [1,0.975,0.95,0.932],
"y": [4.2,4.3,4.412,4.5]
}]