Version 2013-01-06 00:04:43
Inspired by a question from Diego Pujoni on r-sig-mixed, and (a few days later, presumably because there was no response) on r-help.
The parameters of the problem are that the intercept (0-day value) is known to be identical across groups and individuals.
library(nlme)
library(lme4)
library(R2admb)
## library(reshape2)
library(plyr)
library(ggplot2)
theme_set(theme_bw())
library(bbmle) ## for AICtab
library(splines) ## for profiles
Data:
d <- data.frame(Day = rep(c(0,2,4,6,8,10),each=9),
Group = rep(c("c","c","c","t1","t1","t1","t2","t2","t2"),6),
Individual = rep(1:9,6),
X = c(0.71,0.72,0.71,0.72,0.72,0.72,0.70,0.69,0.70,0.72,0.72,
0.71,0.72,0.72,0.71,0.71,0.70,0.71,0.73,0.73,0.69,0.74,
0.69,0.73,0.67,0.71,0.69,0.71,0.71,0.72,0.70,0.71,0.70,
0.52,0.64,0.60,0.70,0.73,0.73,0.67,0.66,0.71,0.47,0.56,
0.54,0.65,0.73,0.73,0.67,0.71,0.58,0.44,0.52,0.58))
Just for the heck of it, plot the data both with lattice and with ggplot2.
library(lattice)
xyplot(jitter(X)~Day, groups=Group, data=d,type=c("a","p"),
auto.key=list(space="right"))
ggplot version has two small advantages:
1. Add lines both by
individual and group average [should be as easy with stat_summary as with
type="a" in xyplot, but there is a bug in the latest ggplot
version …]);
2. Adjust point size rather than jittering to
visualize overlapping points.
(Both of these would of course be possible with a custom panel.xyplot …)
## have to aggregate by hand: stat_summary bug
d2 <- ddply(d,c("Day","Group"),summarise,X=mean(X))
(g1 <- ggplot(d,aes(x=Day,y=X,colour=Group))+
stat_sum(aes(size=factor(..n..)),alpha=0.5)+
scale_size_discrete(range=c(2,5),name="npts")+
geom_line(aes(group=Individual),alpha=0.5)+
geom_line(data=d2,lwd=2,alpha=0.8)+
scale_colour_brewer(palette="Dark2"))
The main conclusions from these pictures are that (1) we probably ought
to be using a nonlinear rather than a linear model; (2) there might be
some heteroscedasticity (larger variance at lower means, as though there
is a “ceiling” to the data at \( \approx X=0.7 \)); it does look as though
there could be among-individual variation (based especially on the t2
data, where the individual curves are approximately parallel). However,
we'll also try linear fits for illustration (even though they won't be
very good):
nlmeLinear fits with lme fail:
LME <- lme(X ~ 1, random = ~Day|Individual, data=d)
## Error: nlminb problem, convergence error code = 1 message = iteration
## limit reached without convergence (10)
If we run this with control=lmeControl(msVerbose=TRUE))) we get a lot
of output, ending with:
47: -65.306481: 5.38940 0.705107 179.050
48: -65.306489: 5.42212 0.705107 184.783
49: -65.306493: 5.45375 0.705106 190.516
50: -65.306499: 5.47352 0.705104 194.382
50: -65.306499: 5.47352 0.705104 194.382
Unsurprisingly, a more complex
model allowing for Group*Day effects fails too:
LME1 <- lme(X ~ Group*Day, random = ~Day|Individual, data=d)
## Error: nlminb problem, convergence error code = 1 message = iteration
## limit reached without convergence (10)
I tried to fit a non-linear model using SSfpl, a self-starting
four-parameter logistic model (with parameters left-asymptote,
right-asymptote, midpoint, scale parameter). This works fine
for an nls fit, giving reasonable results:
nlsfit1 <- nls(X ~ SSfpl(Day, asymp.L, asymp.R, xmid, scale),data=d)
coef(nlsfit1)
## asymp.L asymp.R xmid scale
## 0.7138 0.6248 6.0999 0.9863
Can use gnls to fit group-level differences (although for some reason I
need to specify starting values, even though the help file would lead me
to believe I don't have to … perhaps I do when params is specified?)
My first attempt is apparently a little too ambitious for gnls:
svec <- list(asymp.L=0.7,
asymp.R=c(0.6,0,0),
xmid=c(5,0,0),
scale=c(1,0,0))
gnlsfit1 <- gnls(
X ~ SSfpl(Day, asymp.L, asymp.R, xmid, scale),
data=d,
params=list(asymp.R+scale+xmid~Group,asymp.L~1),
start=svec)
## Error: step halving factor reduced below minimum in NLS step
But I can succeed if I allow only asymp.R to vary among groups:
svec2 <- list(asymp.L=0.7,
asymp.R=c(0.6,0,0),
xmid=6,
scale=1)
gnlsfit2 <- gnls(X ~ SSfpl(Day, asymp.L, asymp.R, xmid, scale),data=d,
params=list(asymp.R~Group,asymp.L++scale+xmid~1),
start=svec2)
Plot predicted values:
predframe <- with(d,expand.grid(Day=unique(Day),Group=levels(Group)))
predframe$X <- predict(gnlsfit2,newdata=predframe)
g1 + geom_line(data=predframe,linetype=3,lwd=2,alpha=0.8)
These look pretty good (it would be nice to get confidence intervals too, but that's a little bit too much of a pain for right now – need to use either delta method or bootstrapping).
dp <- data.frame(d,res=resid(gnlsfit2),fitted=fitted(gnlsfit2))
(diagplot1 <- ggplot(dp,aes(x=factor(Individual),
y=res,colour=Group))+
geom_boxplot(outlier.colour=NULL)+
scale_colour_brewer(palette="Dark2"))
With the exception of individual #7 there's not a lot of evidence for among-individual variation … if we wanted an excuse to ignore among-individual variation we could use
anova(lm(res~Individual,data=dp))
## Analysis of Variance Table
##
## Response: res
## Df Sum Sq Mean Sq F value Pr(>F)
## Individual 1 0.0000 0.000012 0.01 0.91
## Residuals 52 0.0444 0.000854
(whatever philosophical issues this raises about using a large \( p \)-value to accept the null hypothesis that among-individual variation is absent …)
More general diagnostic plot – residual vs. fitted, with points from the same individual connected with lines. There is some hint of decreasing variance with increasing mean.
(diagplot2 <- ggplot(dp,aes(x=fitted,y=res,colour=Group))+geom_point()+
geom_smooth(aes(group=1),colour="black",method="loess")+
geom_path(aes(group=Individual))+
geom_hline(yintercept=0,lty=2))
I can't use nlme with the more ambitious (three parameters varying
by group) model, but I can if I only allow asymp.R to vary:
nlmefit1 <- nlme(model = X ~ SSfpl(Day, asymp.L, asymp.R, xmid, scale),
fixed = list(asymp.R ~ Group, xmid+scale+asymp.L ~1),
random = asymp.R ~ 1 | Individual,
start = list(fixed=with(svec2,c(asymp.R,xmid,scale,asymp.L))),
data=d)
The estimate of the variance in the right-hand asymptote is non-zero (yay):
VarCorr(nlmefit1)
## Individual = pdLogChol(list(asymp.R ~ 1))
## Variance StdDev
## 0.0009292 0.03048
## Residual 0.0004687 0.02165
Adding the random effects doesn't change the parameters much at all:
print(mm1 <- merge(data.frame(cc=coef(gnlsfit2)),
data.frame(cc=fixef(nlmefit1)),by="row.names"),
digits=4)
## Row.names cc.x cc.y
## 1 asymp.L 0.71508 0.71514
## 2 asymp.R.Groupt1 -0.04125 -0.04081
## 3 asymp.R.Groupt2 -0.20162 -0.19986
## 4 asymp.R.(Intercept) 0.71158 0.71162
## 5 scale 0.85397 0.82750
## 6 xmid 5.65084 5.59617
maxdiff <- max(apply(mm1[,-1],1,function(x) abs(diff(x)/mean(x))))
The biggest proportional difference is 3.1% (in the scale
parameter).
nlmefit2 <- update(nlmefit1,fixed = list(asymp.R+xmid+scale+asymp.L ~1),
start = list(fixed=with(svec2,c(asymp.R[1],xmid,scale,asymp.L))))
We can compare the models via AIC or likelihood ratio test
(AICtab from the bbmle package is not essential,
but gives pretty output):
AICtab(nlmefit1,nlmefit2,weights=TRUE)
## dAIC df weight
## nlmefit1 0.0 8 0.999
## nlmefit2 14.2 6 <0.001
anova(nlmefit1,nlmefit2)
## Model df AIC BIC logLik Test L.Ratio p-value
## nlmefit1 1 8 -229.3 -213.4 122.6
## nlmefit2 2 6 -215.1 -203.1 113.5 1 vs 2 18.22 1e-04
It would be nice to do an \( F \) test rather than LRT (i.e. accounting for the finite-size correction), but it's a little bit more work (and probably not really necessary since the answer is so strong).
The
devdiff <- -2*c(logLik(nlmefit2)-logLik(nlmefit1))
pchisq(devdiff,df=2,lower.tail=FALSE)
##
## 0.0001106
## F-test with very large denominator:
pf(devdiff/2,df1=2,df2=1000000,lower.tail=FALSE)
##
## 0.0001106
## Value Std.Error DF t.value Pr(|t|>1)
## asymp.R.(Intercept) 0.7116 0.0206 40.0000 34.47 < 2e-16 ***
## asymp.R.Groupt1 -0.0408 0.0292 40.0000 -1.40 0.16990
## asymp.R.Groupt2 -0.1999 0.0301 40.0000 -6.65 5.8e-08 ***
## xmid 5.5962 0.2382 40.0000 23.49 < 2e-16 ***
## scale 0.8275 0.2175 40.0000 3.81 0.00048 ***
## asymp.L 0.7151 0.0051 40.0000 140.31 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
We don't really know the relevant denominator df, but the summary above suggests the denominator df is 40 (based on the usual type of classical-block-design analogy for df counting, see Pinheiro and Bates 2000 or the glmm wikidot faq).
pf(devdiff/2,df1=2,df2=40,lower.tail=FALSE)
##
## 0.0005493
With help from Reinhold Kliegl (see below) I think I can
now get the correct model specification for nlmer, but
I can't find a convenient syntax for
non-trivial fixed effect modeling
(i.e. in this case allowing some of the
parameters to vary by group) – and when I do construct
the correct syntax, nlmer seems too fragile to get the answer.
This StackOverflow question asks
the same question (syntax, not overcoming fragility).
The basic RE model (without group effects) works fine:
nlmerfit1 <- nlmer(
X ~ SSfpl(Day, asymp.L, asymp.R, xmid, scale) ~
asymp.R|Individual,
start = list(nlpars=c(asymp.L=0.7,
asymp.R=0.6,xmid=5,scale=1)),data=d)
With help from Reinhold Kliegl I got a bit farther (only to
be frustrated again: see below). As I understand it, one
simply has to construct one's own function that encapsulates
the fixed-effect structure; to be used with nlmer, the function
also needs to compute the gradient with respect to the fixed-effect
parameters. This is a bit of a pain, but can be slightly automated
with some hackery by modifying the function generated
by the deriv function:
mm <- model.matrix(~Group,data=d)
grp2 <- mm[,2]
grp3 <- mm[,3]
## or just grp2=as.numeric(d$Group)==2 etc.
nparams <- c("A",paste0("B",0:2),"xmid","scale")
fpl <- deriv(~A+((B0+B1*grp2+B2*grp3)-A)/(1+exp((x-xmid)/scale)),
nparams,
function.arg=c("x",nparams))
attr(fpl,"pnames") <- nparams
tmpstr <- deparse(fpl)
L1 <- grep("^ +\\.value +<-",tmpstr)
L2 <- grep("^ +attr\\(\\.value",tmpstr)
tmpstr2 <- c(tmpstr[1:L1],
paste0(".actualArgs <- as.list(match.call()[",
deparse(nparams),"])"),
tmpstr[(L1+1):(L2-1)],
"dimnames(.grad) <- list(NULL, .actualArgs)",
tmpstr[L2:length(tmpstr)])
fpl <- eval(parse(text=tmpstr2))
Attempt the fit:
nlmerfit2 <- nlmer(
X ~ fpl(Day, asymp.L, asymp.R1, asymp.R2, asymp.R3, xmid, scale) ~
asymp.R1|Individual,
start = list(nlpars=c(asymp.L=0.7,
asymp.R1=0.6,asymp.R2=0,asymp.R3=0,xmid=5,scale=1)),data=d)
## Error: step factor reduced below 0.001 without reducing pwrss
Fails (but I think this is due to the fragility of nlmer rather
than anything fundamentally wrong with the setup). To be sure
I should generate a larger, artificial data set along the same lines
and see if I can get it to work.
Right now we are using the development version of lme4, i.e. version 0.99999911.0.
We can get an answer (but not a sensible answer) with the stable (lme4.0) version.
detach("package:lme4",unload=TRUE)
library("lme4.0")
## Attaching package: 'lme4.0'
## The following object is masked from 'package:nlme':
##
## lmList, VarCorr
packageVersion("lme4.0")
## [1] '0.999999.1'
nlmerfit2 <- nlmer(
X ~ fpl(Day, asymp.L, asymp.R1, asymp.R2, asymp.R3, xmid, scale) ~
asymp.R1|Individual,
start = c(asymp.L=0.7,
asymp.R1=0.6,asymp.R2=0,asymp.R3=0,xmid=5,scale=1),data=d)
The results are silly!
fixef(nlmerfit2)
## asymp.L asymp.R1 asymp.R2 asymp.R3 xmid scale
## -156.16916 45.81031 -0.01931 -0.12157 4428.96516 3548.54685
range(predict(nlmerfit2))
## [1] -156.2 4429.0
I can't be sure (short of getting the Word of God from DB) that
there is no easier way to do fixed effects in nlmer, but I have stared
at the source code and spent several hours struggling with
various permutations to do this in a more “built-in” way, without
success.
We can also use AD Model Builder to tackle this problem. It seems (see below) that it can handle the more complex models, such as fitting group effects of more of the parameters; to my surprise, some of these models turn out to be informative (I expected that we were more or less at the limits of the signal in the data).
Partly due to my lower level of familiarity with ADMB, this was a bit of a slog, and I only succeeded in the end by building up via gradual steps.
First try it without random effects, grouping variables, etc.
(i.e. equivalent to nls fit above).
setup_admb()
## [1] "/usr/local/admb"
admb_model0 <- "
PARAMETER_SECTION
vector pred(1,nobs) // predicted values
PROCEDURE_SECTION
pred = asympL+(asympR-asympL)/(1.0+exp(-(Day-xmid)/scale));
f = regression(X,pred);
"
writeLines(admb_model0,con="algae0.tpl")
## set up data: adjust names, etc.
d0 <- subset(d,select=c(Day,X))
d0 <- c(list(nobs=nrow(d)),as.list(d0))
## starting values: adjust names, add values
svec3 <- svec2
names(svec3) <- gsub("\\.","",names(svec3)) ## remove dots
svec3$asympR <- 0.6 ## single value
## run
admbfit0 <- do_admb("algae0",
data=d0,
params=svec3,
run.opts=run.control(checkparam="write",
checkdata="write"))
Works fine:
## admb nls
## asympL 0.713782 0.713782
## asympR 0.624793 0.624793
## xmid 6.099903 6.099905
## scale 0.986272 0.986270
Now try with fixed effects grouping (equivalent
to gnlsfit2 above), using the dummy variables
constructed above (could also use if statements, or a vector
of R values in a for loop with R[Group[i]],
or (best choice) pass a model matrix for R …). We have
to use elem_div tather than /
for elementwise division of two vectors.
admb_model1 <- "
PARAMETER_SECTION
vector pred(1,nobs) // predicted values
vector Rval(1,nobs) //
PROCEDURE_SECTION
Rval = asympR1+asympR2*grp2+asympR3*grp3;
pred = asympL+elem_div(Rval-asympL,1.0+exp(-(Day-xmid)/scale));
f = regression(X,pred);
"
writeLines(admb_model1,con="algae1.tpl")
d1 <- c(d0,list(grp2=grp2,grp3=grp3))
svec4 <- c(rename(svec3,c(asympR="asympR1")),
list(asympR2=0,asympR3=0))
admbfit1 <- do_admb("algae1",
data=d1,
params=svec4,
run.opts=run.control(checkparam="write",
checkdata="write"))
Trying it with a model matrix instead:
admb_model1B <- "
PARAMETER_SECTION
vector pred(1,nobs) // predicted values
vector Rval(1,nobs) //
PROCEDURE_SECTION
Rval = XR*Rvec;
pred = asympL+elem_div(Rval-asympL,1.0+exp(-(Day-xmid)/scale));
f = regression(X,pred);
"
writeLines(admb_model1B,con="algae1B.tpl")
d1B <- c(d0,list(XR=mm))
svec4B <- rename(svec3,c(asympR="Rvec"))
svec4B$Rvec <- c(0.6,0,0)
admbfit1B <- do_admb("algae1B",
data=d1B,
params=svec4B,
run.opts=run.control(checkparam="write",
checkdata="write"))
Works, with identical parameters (except for order/names), of course.
Now adding the random effects.
The regression function is not fully implemented for random effects
models
(although this should be fixed in an upcoming version), so
we subtitute the formula \( n/2 \log(\text{RSS}/n) \), where RSS
is the residual sum of squares (coded as norm2(X-pred)):
admb_model2 <- "
PARAMETER_SECTION
vector pred(1,nobs) // predicted values
vector Rval(1,nobs) //
PROCEDURE_SECTION
Rval = XR*Rvec+Rsigma*(Z*Ru);
pred = asympL+elem_div(Rval-asympL,1.0+exp(-(Day-xmid)/scale));
f = 0.5*nobs*log(norm2(X-pred)/nobs)+norm2(Ru);
"
writeLines(admb_model2,con="algae2.tpl")
d2 <- c(d0,list(XR=mm,Z=model.matrix(~factor(Individual)-1,data=d)))
svec5 <- c(svec4B,list(Rsigma=list(0.1,bounds=c(0,10))))
admbfit2 <- do_admb("algae2",
data=d2,
params=svec5,
re=list(Ru=max(d$Individual)),
run.opts=run.control(checkparam="write",
checkdata="write"))
## [1] TRUE
Since ADMB doesn't handle sparse matrices, and doesn't penalize
looping, it would be slightly more efficient to implement the
random effects as for (i=1; i<=nobs; i++) Rval[i] += Rsigma*Ru[Group[i]],
but I'm being lazy/I like the compactness and extensibility of
the matrix representation …
Now we can finally test for fixed-effect differences in parameters
other than R:
admb_model3 <- "
PARAMETER_SECTION
vector pred(1,nobs) // predicted values
vector Rval(1,nobs) //
vector scaleval(1,nobs)
vector xmidval(1,nobs)
sdreport_vector opred(1,nobs) // output predictions
PROCEDURE_SECTION
Rval = XR*Rvec+Rsigma*(Z*Ru);
scaleval = Xscale*scalevec;
xmidval = Xxmid*xmidvec;
pred = asympL+elem_div(Rval-asympL,
1.0+exp(-elem_div(Day-xmidval,scaleval)));
opred=pred;
f = 0.5*nobs*log(norm2(X-pred)/nobs)+norm2(Ru);
"
writeLines(admb_model3,con="algae3.tpl")
d3 <- c(d2,list(Xscale=mm,Xxmid=mm))
svec6 <- rename(svec5,c(xmid="xmidvec",scale="scalevec"))
svec6$xmidvec <- c(svec6$xmidvec,0,0)
svec6$scalevec <- c(svec6$scalevec,0,0)
admbfit3 <- do_admb("algae3",
data=d3,
params=svec6,
re=list(Ru=max(d$Individual)),
run.opts=run.control(checkparam="write",
checkdata="write"))
Results:
AICtab(admbfit0,admbfit1,admbfit2,admbfit3,weights=TRUE)
## dAIC df weight
## admbfit3 0.0 11 0.99761
## admbfit2 12.1 7 0.00239
## admbfit1 34.9 6 < 0.001
## admbfit0 111.3 4 < 0.001
summary(admbfit3)
## Model file: algae3_gen
## Negative log-likelihood: -214.2 AIC: -450.5
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## asympL 0.71222 0.00341 208.62 < 2e-16 ***
## Rvec.1 0.70342 0.03399 20.69 < 2e-16 ***
## Rvec.2 -0.06823 0.06525 -1.05 0.2957
## Rvec.3 -0.18892 0.04592 -4.11 3.9e-05 ***
## xmidvec.1 9.31528 1.59892 5.83 5.7e-09 ***
## xmidvec.2 -0.27981 1.56654 -0.18 0.8582
## xmidvec.3 -3.79499 1.60802 -2.36 0.0183 *
## scalevec.1 0.96047 0.70150 1.37 0.1709
## scalevec.2 -0.31566 0.87108 -0.36 0.7171
## scalevec.3 -0.15714 0.71820 -0.22 0.8268
## Rsigma 0.07297 0.02290 3.19 0.0014 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
There's a very large AIC difference. Likelihood ratio/\( F \) tests
as shown above for the nlme fits are left as an exercise …
For the plot, it's nice to redo the fit with the parameters specified per-group, rather than as baseline+contrast:
d4 <- d3
d4$XR <- d4$Xscale <- d4$Xxmid <- model.matrix(~Group-1,data=d)
svec7 <- svec6
for (i in c("Rvec","xmidvec","scalevec"))
svec7[[i]] <- rep (svec7[[i]][1],3)
admbfit3B <- do_admb("algae3",
data=d4,
params=svec7,
re=list(Ru=max(d$Individual)),
run.opts=run.control(checkparam="write",
checkdata="write"
))
library(coefplot2)
## n.b. need latest version (1.3.3.3)
cc <- coeftab(admbfit3B)[1:11,] ## omit RE parameters
coefplot2(list(cc),intercept=TRUE)
Now we're annoyed about the scaling, so (with some fussing) we can replot the group-varying parameters in separate panels:
Diagnostic plots:
## drop conditional mode/individual-R estimates
admb3pred <- coef(admbfit3,"extra")[-(1:9)]
admb3res <- admb3pred-d$X
dp2 <- data.frame(d,res=admb3res,fitted=admb3pred)
diagplot1 %+% dp2
Perhaps a hint of greater variance in the two treatment groups?
Fitted vs. residual:
diagplot2 %+% dp2
Overlay predictions (dashed lines):
g1 + geom_line(data=dp2,aes(y=fitted,group=Individual),lty=2)
It would be nice but a bit more tedious to generate smooth
prediction curves (i.e. for intermediate Day values).
xmidvec is smaller) as well as dropping farther
(Rvec is lower).Computing a likelihood profile for \( \sigma^2_R \) is not too hard to code (now that I've implemented phase support in R2admb, which allows us to fix the values of selected parameters by making their phases negative), but is a bit of a pain to run: it's slow, and the computation craps out with non-positive-definite matrices at several points right around the lower confidence bound; I ran another set of values to try to cover the region adequately (grubby details hidden).
svec8 <- svec6
Rsigmavec <- seq(0,0.24,by=0.01)
fitfun <- function(Rsigma) {
## cat(".")
svec8$Rsigma <- list(Rsigma,phase=-1)
try(do_admb("algae3",
data=d3,
params=svec8,
re=list(Ru=max(d$Individual)),
run.opts=run.control(checkparam="write",
checkdata="write")),silent=TRUE)
}
admb_Rsigma_prof <- lapply(Rsigmavec,fitfun)
## try to fill in holes
Rsigmavec2 <- seq(0.032,0.048,by=0.002)
admb_Rsigma_prof2 <- lapply(Rsigmavec2,fitfun)
The profiled points with the interpolation spline and the likelihood ratio test cutoff line:
The 95% profile confidence intervals on \( \sigma^2_R \) (derived from inverting an interpolation spline fitted to these points) are {0.0386,0.2169}.
I haven't computed it, but the transformed profile plot (on a \( y \) scale corresponding to the square-root deviation of the deviance from the minimum deviance, so that a quadratic profile would be a symmetric V) shows that the quadratic approximation would be fairly bad for this case …
ggplot(na.omit(proftab12),aes(Rsigma,sqrt(2*(NLL-min(NLL,na.rm=TRUE)))))+
geom_point()+labs(x=~{sigma^2}[R],y=~sqrt(Delta*"deviance"))+
geom_hline(yintercept=1.96,colour="red")+
geom_line(data=sdat,colour="blue",alpha=0.7)
xmid and scale
parameters as well, although I think that would be a little ambitious
for this data set.X.