####################
# Code for variogram
####################

# fit lm1, a linear model with the respective structure

r <- resid(lm1) # get residuals from linear model

t <- rats$logT # get time points
id <- rats$SUBJECT # get IDs


half.squared<- outer(r, r, FUN = function(r1, r2) (r1 - r2)^2/2) # half squared residuals for all observations
distances <- outer(t, t, FUN = function(t1, t2) abs(t1 - t2)) # absolute distances for all observations
samei <- outer(id, id, FUN = function(i1, i2) (i1 == i2)) # get index matrix specifying if same subject
nondiag <- 1 - diag(length(r)) # matrix with 0 on the diagonal and 1 for non-diagonal values as 0 on diagonal
u <- as.vector(distances[samei & nondiag]) 
v <- as.vector(half.squared[samei & nondiag])
totvar <- mean(half.squared[!samei]) # total variance = \sigma^2 + d^2+ \tau^2

# smooth and plot
smooth_mean <- lowess(u,v,iter=0) # note: iter = 0 such that not the median is approximated, see ?lowess
par(mar=c(5.1, 5.1, 4.1, 2.1))
plot(u, v, pch = 16, cex = 0.5, col = 8,ylim=c(0,totvar),xlim=c(0,2.005),
     xlab="distance of logT (u)",ylab=expression(paste(hat("v"),"(u)")),main="Empirical Semi-Variogram")
lines(smooth_mean, lwd = 2,t="b")
abline(h = totvar, col = 2, lwd = 2, lty = 2)