Solutions

1

Consider the ToothGrowth data.

1. Compute using apply, lapply etc the average growth (variable len) for the two different supplements (variable supp).
2. Compute using apply, lapply etc the average growth (variable len) for each combination of supplement (variable supp) and dose (variable dose).

attach(ToothGrowth)
avg_s = sapply(c("OJ", "VC"), function(s) {mean(len[supp == s])})
avg_s

OJ

20.6633333333333

VC

16.9633333333333

avg_sd = sapply(c(0.5, 1.0, 2.0), function(d) {sapply(c("OJ", "VC"), function(s) {mean(len[supp == s & dose == d])})})
avg_sd
# OR
avg_sd = mapply(function(s, d) {mean(len[supp == s & dose == d])}, 
                 rep(c("OJ", "VC"), 3), rep(c(0.5, 1.0, 2.0), each = 2))
avg_sd

OJ	13.23	22.70	26.06
VC	7.98	16.77	26.14

OJ

13.23

VC

7.98

OJ

22.7

VC

16.77

OJ

26.06

VC

26.14

2

Consider the topo data provided by the package MASS containing the altitude z at 52 coordinates x and y.

1. Use for loops efficiently to compute the matrix of pairwise distances between the locations given by the x, y coordinates. Compute only the lower triangular part of the matrix.
2. The spherical correlation structure is given by the function $$r(d; \theta) = \begin{cases} 1 - 1.5 \dfrac{d}{\theta} + 0.5 \biggl( \dfrac{d}{\theta} \biggr)^3 & \text{if $d < \theta$,} \\ 0 & \text{otherwise.} \end{cases}$$ where $d$ is the distance between two coordinates and $\theta$ is a parameter that controls the strength of the correlation. Write efficient parallel code to compute the correlation matrix for the sampling coordinates of the topo data set at $\theta = 3$.

library(MASS)
attach(topo)
n = nrow(topo)
distmat = matrix(0, n, n) # Initialise memory for the matrix of distances
for (j in 1:(n-1)) { # Loop across columns
  for (i in (j+1):n) { # Loop across the lower triangular part
    distmat[i, j] = sqrt((x[i] - x[j])^2 + (y[i] - y[j])^2)
  }
}
corrmat = matrix(0, n, n) # Initialise memory for the correlation matrix
theta = 3
rowindx = matrix(1:n, n, n)
lowtri = rowindx > t(rowindx) # Identify the lower triangular part of the matrix
idisttheta = distmat < theta # Identify elements < theta
disttheta = distmat[idisttheta & lowtri]/theta
corrmat[idisttheta & lowtri] = 1 - 1.5*disttheta + .5*disttheta^3

The following objects are masked _by_ .GlobalEnv:
    x, z