# This file contains code to generate the model estimates in
# An Empirical Model of Slope Ratio Comparisons
# Justin Talbot, John Gerth, and Pat Hanrahan
# InfoVis 2012

n <- read.csv("results_study1.csv", header=TRUE)

n$midangle <- ((n$theta1+n$theta2)/2)/pi*180

n$slope.ratio <- n$actual
n$angle.ratio <- (n$theta2/n$theta1)*100

Height <- n[n$strategy=="Height",]
Angle  <- n[n$strategy=="Angle",]

### Height

# Compute the height model (parameter gamma in 6.2)
# offset() forces slope.ratio to have a coefficient of 1, since
# our model (Equation 2) doesn't have a coefficient on this term
height <- lm(response~offset(slope.ratio)+1, data=Height)
summary(height)


### Angle

# Compute the angle model (paramters mu and beta in 6.2)
# Note that the angle model uses  the true angle ratio, not the true slope ratio.
# Again, our model doesn't have a coefficient on the true angle ratio term,
# so we use offset()
angle <- lm(response~offset(angle.ratio)+midangle+1, data=Angle)
summary(angle)

# With an offset in the model formula, R doesn't compute the R^2 I want.
# I want to compare my model to the model response~1. 
# So, I'll compute it manually.
tss <- sum((Angle$response-mean(Angle$response))^2)
rss <- sum(summary(angle)$residuals^2)
1-rss/tss

# Try to fit models with only one of the two components
# Can't look at R^2 for this first model since it
# doesn't have an intercept, but anova will work fine.
# Can't use anova for the second model, since it doesn't
# like the missing offset, so just look at the rss.
angle.nomidangle <- lm(response~offset(angle.ratio)-1, data=Angle)
angle.noangleratio <- lm(response~midangle+1, data=Angle)

anova(angle.nomidangle, angle)
anova(angle.noangleratio, angle)

# Try to fit model with a coefficient on the true angle
angle.coef <- lm(response~angle.ratio+midangle+1, data=Angle)
summary(angle.coef)
anova(angle, angle.coef)

# Test presence of interaction
# (slightly significant, but magnitude is too small to be worthwhile)
angle.interaction <- lm(response~angle.ratio*midangle+1, data=Angle)
summary(angle.interaction)
anova(angle, angle.interaction)

