# square root models baseline <- read.csv ("IPM_BASELINE_R2_032006.csv") attach.all (baseline) sqrt.roach1 <- sqrt (roach1) sqrt.roach2 <- sqrt (roach2) M1 <- lm (sqrt.roach2 ~ sqrt.roach1 + treatment) display (M1) M1a <- lm (sqrt.roach2 ~ sqrt.roach1 + treatment + factor(building)) display (M1a) M2 <- lmer (sqrt.roach2 ~ sqrt.roach1 + treatment + (1 | building)) display (M2) M3 <- lmer (sqrt.roach2 ~ sqrt.roach1 + treatment + senior + (1 | building)) display (M3) # effects of -1.27 # if it were 10 under the control print ((sqrt(10) - 1.27)^2) # start at 50 print ((sqrt(50) - 1.27)^2) # start at 100 print ((sqrt(100) - 1.27)^2) # start at 1000 print ((sqrt(1000) - 1.27)^2) # effects of -1.27 - 2*.59 = -2.45 # if it were 10 under the control print ((sqrt(10) - 2.45)^2) # etc sqrt.diff <- function (start, change){ return (new.value=(sqrt(start) + change)^2) } diffs <- function (beta.hat){ cat ("assuming beta.hat =", beta.hat, "\n") for (start in c(0,5,10,50,100)){ cat ("instead of", start, ", it's", sqrt.diff (start, beta.hat), "(on average)\n") } } diffs (fixef(M3)["treatment"]) diffs (fixef(M3)["treatment"] + 2*se.fixef(M3)["treatment"]) diffs (fixef(M3)["treatment"] - 2*se.fixef(M3)["treatment"]) # redo in terms of roach1 b.hat <- fixef (M3) s.hat <- sqrt (sigma.hat(M3)\$sigma\$data^2 + sigma.hat(M3)\$sigma\$building^2) diffs2 <- function (senior.value, beta, sigma){ r1 <- median(roach1[senior==senior.value])*c(.5,1,2) N <- length (r1) pred.T <- rep (NA, N) pred.C <- rep (NA, N) for (i in 1:N){ x.T <- cbind (1, sqrt(r1[i]), 1, senior.value) # treatment x.C <- cbind (1, sqrt(r1[i]), 0, senior.value) # control pred.T[i] <- mean (rnorm (500, x.T %*% beta, sigma)^2) pred.C[i] <- mean (rnorm (500, x.C %*% beta, sigma)^2) } output <- cbind (r1, pred.T, pred.C, pred.C - pred.T) dimnames(output) <- list (c ("low", "median", "high"), c ("r1", "pred.T", "pred.C", "pred.T - pred.C")) return (output) } diffs2 (0, b.hat, s.hat) diffs2 (1, b.hat, s.hat) # now embed in a simulation n.sims <- 200 beta.sim <- sim (M3, n.sims)\$unmodeled diffs2.sim <- function (senior.value, beta.sim, sigma){ n.sims <- nrow (beta.sim) output <- array (NA, c(n.sims,3,4)) for (s in 1:n.sims){ output[s,,] <- diffs2 (senior.value, beta.sim[s,], sigma) } for (j in 1:3){ r1 <- mean (output[,j,1]) mean.T <- mean (output[,j,2]) mean.C <- mean (output[,j,3]) mean.diff <- mean (output[,j,4]) i50.diff <- quantile (output[,j,4], c(.25, .75)) i95.diff <- quantile (output[,j,4], c(.025, .975)) cat ("Senior.value = ", senior.value, ": If roach1 = ", fround(r1,1), ", then mean.T = ", fround(mean.T,1), ", mean.C = ", fround(mean.C,1), ",\n and 50% and 95% conf intervals for diff are [", fround(i50.diff[1], 1), ", ", fround(i50.diff[2], 1), "] and [", fround(i95.diff[1], 1), ", ", fround(i95.diff[2], 1), "]\n", sep="") } } diffs2.sim (0, beta.sim, s.hat) diffs2.sim (1, beta.sim, s.hat)