### ICC : Example 1

dat1 <-read.table("lecture12ex1.csv", sep=",", header=TRUE)
head(dat1)
psych::describe(dat1)

with(dat1, cor(menneuro, womenneuro))
with(dat1, cor(mendistress, womendistress))

# ICC : Example 2

dat2 <-read.table("lecture12ex2.csv", sep=",", header=TRUE)
head(dat2)
psych::describe(dat2)

library(lme4)
out2a <- lmer(extro ~ 1 + (1|coupleid), data=dat2, REML=FALSE)
summary(out2a) 
out2b <- lmer(neuro ~ 1 + (1|coupleid), data=dat2, REML=FALSE)
summary(out2b) 

library(nlme)
out2c <- gls(extro ~ 1, data = dat2, 
             correlation = corCompSymm(form = ~1|coupleid))
summary(out2c)
out2d <- gls(neuro ~ 1, data = dat2, 
             correlation = corCompSymm(form = ~1|coupleid))
summary(out2d)

result <- rep(NA, 8)
dimensions <- 3:10
for(i in seq_along(dimensions)) {
  k <- dimensions[i]
  f <- function(r) {
    R <- matrix(r, k, k)
    diag(R) <- 1
    solve(eigen(R)$values[k])
  }
  result[i] <- uniroot(f, lower=-0.99, upper=-0.01)$root
}
plot(dimensions, result, type='l', ylim=c(-0.6, 0.001), 
     xlab="Number of Members in Groups", 
     ylab="Maximum Possible Negative ICC")
points(dimensions, result)

dat2a <- dat2[order(dat2$coupleid, dat2$subjectid),]
head(dat2a)
dat2b <- dat2[order(dat2$coupleid, -dat2$subjectid),]
head(dat2b)
dat2a$extropartner <- dat2b$extro
dat2a$neuropartner <- dat2b$neuro
head(dat2a)
with(dat2a, cor(extro, extropartner))
with(dat2a, cor(neuro, neuropartner))

# Distinguishabiltiy

library(lavaan)
dat1 <-read.table("lecture12ex1.csv", sep=",", header=TRUE)
scriptparent <- '
# Mean
mendistress ~ 1
womendistress ~ 1
menneuro ~ 1
womenneuro ~ 1

# Variance
mendistress ~~ mendistress
womendistress ~~ womendistress
menneuro ~~ menneuro
womenneuro ~~ womenneuro

# Covariances
mendistress ~~ womendistress 
menneuro ~~ womenneuro
mendistress ~~ menneuro 
womendistress ~~ womenneuro
mendistress ~~ womenneuro
womendistress ~~ menneuro
'

scriptnested <- '
# Mean
mendistress ~ m1*1
womendistress ~ m1*1
menneuro ~ m2*1
womenneuro ~ m2*1

# Variance
mendistress ~~ v1*mendistress
womendistress ~~ v1*womendistress
menneuro ~~ v2*menneuro
womenneuro ~~ v2*womenneuro

# Covariances
mendistress ~~ dis*womendistress 
menneuro ~~ neuro*womenneuro
mendistress ~~ same*menneuro 
womendistress ~~ same*womenneuro
mendistress ~~ diff*womenneuro
womendistress ~~ diff*menneuro 
'

outparent <- sem(scriptparent, data=dat1)
outnested <- sem(scriptnested, data=dat1)
anova(outnested, outparent)
summary(outparent, standardize=TRUE)

# APIM : Distinguishability

dat1 <-read.table("lecture12ex1.csv", sep=",", header=TRUE)

library(lavaan)
script1 <- '
mendistress ~ menneuro + womenneuro + lengthrela
womendistress ~ menneuro + womenneuro + lengthrela
mendistress ~~ womendistress
'
out1 <- sem(script1, data=dat1)
summary(out1, standardize=TRUE, rsquare=TRUE)

script2 <- '
mendistress ~ am*menneuro + pmw*womenneuro + cm*lengthrela
womendistress ~ pwm*menneuro + aw*womenneuro + cw*lengthrela
mendistress ~~ womendistress

# Define New Parameters
apmen := am - pmw                         # Compare actor and partner on mendistress
apwomen := aw - pwm                       # Compare actor and partner on womendistress
diffa := am - aw                          # Compare actor effects on each role
diffp := pmw - pwm                        # Compare partner effects on each role
diffmenneuro := am - pwm                  # Compare the effects on menneuro on men and women distress
diffwomenneuro := pmw - aw                # Compare the effects on womenneuro on men and women distress
difflengthrela := cm - cw                 # Compare the effects on lengthrela on men and women distress
aveactor := (am + aw) / 2                 # Average actor effects
avepartner := (pmw + pwm) / 2             # Average partner effects
diffactorpartner := aveactor - avepartner # Compare actor and partner effects
'
out2 <- sem(script2, data=dat1)
summary(out2)

# APIM : Indistinguishable Dyads

dat3 <- read.table("lecture12ex3.csv", sep=",", header=TRUE)
psych::describe(dat3)

dat3a <- dat3[order(dat3$coupleid, dat3$subjectid),]
dat3b <- dat3[order(dat3$coupleid, -dat3$subjectid),]
dat3a$housepartner <- dat3b$house
dat3a$avehouse <- ave(dat3a$house, dat3a$coupleid)

library(nlme)
out3a <- gls(sat ~ 1 + house + housepartner,
            data = dat3a, correlation = corCompSymm(form = ~1|coupleid))
summary(out3a)

sumout3a <- summary(out3a)
psi <- coef(sumout3a)[,1]
vpsi <- vcov(sumout3a)
cvec1 <- c(0, 1, -1) 
cmat <- cbind(cvec1)
est <- t(cmat) %*% psi
se <- sqrt(diag(t(cmat) %*% vpsi %*% cmat))
zval <- est/se
pval <- round(pnorm(-abs(zval)) * 2, 6)
cbind(est, se, zval, pval)

out3b <- gls(sat ~ 1 + house*housepartner,
             data = dat3a, correlation = corCompSymm(form = ~1|coupleid))
summary(out3b)

mean(dat3$house)
sd(dat3$house)

plot(dat3a$house, dat3a$sat, type="n", xlab="Self Housework Hours per Month", 
     ylab="Relationship Satisfaction", main="Actor Effect")
abline(a = 53.56, b = -3.19, col="red")
abline(a = 57.13, b = -2.07, col="green")
abline(a = 60.70, b = -0.95, col="blue")
legend(0.1, 40, c("Low (2.25)", "Medium (3.25)", "High (4.25)"), 
       col=c("red", "green", "blue"), lty=1, title="Partner Housework Hours per Month")

sumout3b <- summary(out3b)
psi <- coef(sumout3b)[,1]
vpsi <- vcov(sumout3b)
cvec1 <- c(1, 0, 2.25, 0) 
cvec2 <- c(1, 0, 3.25, 0) 
cvec3 <- c(1, 0, 4.25, 0) 
cvec4 <- c(0, 1, 0, 2.25) 
cvec5 <- c(0, 1, 0, 3.25) 
cvec6 <- c(0, 1, 0, 4.25) 
cmat <- cbind(cvec1, cvec2, cvec3, cvec4, cvec5, cvec6)
est <- t(cmat) %*% psi
se <- sqrt(diag(t(cmat) %*% vpsi %*% cmat))
zval <- est/se
pval <- round(pnorm(-abs(zval)) * 2, 6)
cbind(est, se, zval, pval)

plot(dat3a$house, dat3a$sat, type="n", xlab="Partner Housework Hours per Month", 
     ylab="Relationship Satisfaction", main="Partner Effect")
abline(a = 32.68, b = 6.09, col="red")
abline(a = 26.97, b = 7.21, col="green")
abline(a = 21.26, b = 8.33, col="blue")
legend(0.1, 90, c("Low (2.25)", "Medium (3.25)", "High (4.25)"), 
       col=c("red", "green", "blue"), lty=1, title="Self Housework Hours per Month")

sumout3b <- summary(out3b)
psi <- coef(sumout3b)[,1]
vpsi <- vcov(sumout3b)
cvec1 <- c(1, 2.25, 0, 0) 
cvec2 <- c(1, 3.25, 0, 0) 
cvec3 <- c(1, 4.25, 0, 0) 
cvec4 <- c(0, 0, 1, 2.25) 
cvec5 <- c(0, 0, 1, 3.25) 
cvec6 <- c(0, 0, 1, 4.25) 
cmat <- cbind(cvec1, cvec2, cvec3, cvec4, cvec5, cvec6)
est <- t(cmat) %*% psi
se <- sqrt(diag(t(cmat) %*% vpsi %*% cmat))
zval <- est/se
pval <- round(pnorm(-abs(zval)) * 2, 6)
cbind(est, se, zval, pval)

out3b <- gls(sat ~ 1 + house*housepartner,
             data = dat3a, correlation = corCompSymm(form = ~1|coupleid))
summary(out3b)

# Rsquare

cor(predict(out3b), dat3$sat)^2

out3c <- lme(sat ~ 1 + house*housepartner, random = ~1 |coupleid, data = dat3a)
library(r2mlm)
r2mlm(out3c)

# Bias 

dat3c <- dat3[!duplicated(dat3$coupleid),]
summary(lm(sat ~ house, data=dat3c))

summary(out3a)