[R] Network meta-analysis, varConstPower in nlme

[Christian Gold]

Dear Thomas Lumley, and R-help list members,

I have read your article "Network meta-analysis for indirect treatment
comparisons" (Statist Med, 2002) with great interest. I found it very
helpful that you included the R code to replicate your analysis;
however, I have had a problem replicating your example and wondered if
you are able to give me a hint. When I use the code from the article:

lme1 <- lme(Y1 ~ trt.B + trt.C + trt.D + trt.E, random = ~ 1 | trtpair,
data=lumley1, var = varConstPower(form=~sigma, fixed=list(power=1)))

I get an error message:

Error in lme(Y1 ~ trt.B + trt.C + trt.D + trt.E, random = ~1 | trtpair,  :
   unused argument(s) (var = list(const = numeric(0), power = numeric(0)))

The problem seems to be in the varConstPower component, but I don't
understand exactly what is wrong. (I have pasted the complete code to
reproduce the result at the end of this email.)
I use R 2.6.1 and nlme 3.1-86.

Any help would be greatly appreciated!

Best,

Christian Gold

-------

lumley1 <- scan()
2 1 -2.428 -1.451
2 5 3.233 3.819
3 2 -0.493 -0.405
1 2 0.577 1.868
3 1 -2.06 -0.709
3 4 1.976 1.196
2 3 1.084 2.055
2 3 1.255 1.278
3 4 0.807 1.058
3 4 0.254 1.672
1 4 3.008 3.634
1 5 3.817 4.063
2 1 -0.755 -0.405
1 3 1.339 2.100

lumley1 <- data.frame(matrix(lumley1, byrow=T, ncol=4))
names(lumley1) <- c("Active", "Control", "Y1", "Y2")

trt <- matrix(0, ncol=5, nrow=nrow(lumley1)); colnames(trt) <-
paste("trt", LETTERS[1:5], sep=".")
for (i in 1:nrow(trt))    {trt[i, lumley1$Active[i]] = 1; trt[i,
lumley1$Control[i]] = -1}
lumley1 <- data.frame(data.frame(trt), lumley1)
lumley1 <- transform(lumley1, trtpair = paste(Active, Control))
sigma <- rep(sqrt(.5), nrow=lumley1)
library(nlme)
lme1 <- lme(Y1 ~ trt.B + trt.C + trt.D + trt.E, random = ~ 1 | trtpair,
data=lumley1, var = varConstPower(form=~sigma, fixed=list(power=1)))
summary(lme1)
lme2 <- lme(Y2 ~ trt.B + trt.C + trt.D + trt.E, random = ~ 1 | trtpair,
data=lumley1), var = varConstPower(form=~sigma, fixed=list(power=1)))
summary(lme2)