Literate R report
2018 February 22
R preliminaries
Before starting the analysis, we ensure the preliminaries: packages; external files; and functions defined ad hoc.
Required files
We ensure that the required files are in the working directory.
"fromajauketal2013.csv", a CSV file of data from Jauk et al 2013."algsvdm.cpp", a C++ source file implementing algorithms from Van Dorp & Mazzuchi 2002.
# list files in working directory
list.files(getwd())## [1] "algsvdm.cpp" "fromjauketal2013.csv" "literaterreport.Rmd"Required packages
We ensure that the required packages are installed in R. Required packages are the following:
Rcpp, for interface with C++.MASS, for fitting data to a marginal beta distribution.segmented, for segmented regression.devtools, for displaying R session info.
# load required packages
library("Rcpp")
library("MASS")
library("segmented")
library("devtools")Functions defined ad hoc
We define a suite of functions for use in our analysis.
The algorithms in Van Dorp & Mazzuchi 2002 are computationally heavy, so we implement them in C++ to optimize performance; otherwise, the duration of the procedure would be prohibitively long. Importing the C++ source file provides function nesurb, our workhorse for matching any pair of quantiles to the corresponding pair of beta shape parameters. See the separate C++ source file for more details.
# use algorithms from Van Dorp & Mazzuchi 2002
Rcpp::sourceCpp("./algsvdm.cpp")We define function nesucheck, which checks the validity of inputs with respect to the model proposed. For instance, the model assumes that \(0 \leq Y \leq 1\); otherwise, carrying out the procedure is erroneous.
# function: ensure that inputs make sense in the model
nesucheck = function(theta,x,y) {
# on theta
if(!is.null(theta)) {
stopifnot(is.numeric(theta)) # numeric
stopifnot(length(theta)==4) # 4 parameters
stopifnot(theta[1]<theta[2]) # floor less than ceiling
stopifnot(theta[3]<theta[4]) # floor less than ceiling
stopifnot(0<theta&theta<1) # between 0 and 1
}
# on y
if(!is.null(y)) {
stopifnot(is.numeric(y)) # numeric
stopifnot(length(x)==length(y)) # same number of observations
stopifnot(0<=y&y<=1) # between 0 and 1
}
# on x
if(!is.null(x)) {
stopifnot(is.numeric(x)) # numeric
}
# valid
return(TRUE)
}We define function nesugetbeta, which obtains the beta shape parameters from \({\theta}\) and \(\boldsymbol{x}\). Uses function nesurb under the hood. Useful for simulation and computing the likelihood.
# function: get beta shape parameters
nesugetbeta = function(theta,x) {
# check compliance to model
nesucheck(theta=theta,x=x,y=NULL)
theta = setNames(as.list(theta),c("a0","b0","a1","b1")) # assign names
# get quantile rank
u = pnorm(x,mean=100,sd=15)
# get quantile
au = with(theta,a0+(a1-a0)*u)
bu = with(theta,b0+(b1-b0)*u)
# get beta shape parameters
ans = nesurb(au,bu) # uses V&M 2002
return(ans)
}We define function nesusimulate, which simulates \({Y \rvert \{X = x\}}\). Set argument keepx to TRUE to have the return object rememeber the predictor x; otherwise the return object is just the simulated response.
# function: simulate response
nesusimulate = function(theta,x,keepx=FALSE) {
# check compliance to model
nesucheck(theta=theta,x=x,y=NULL)
# get beta shape parameters
betashape = nesugetbeta(theta,x)
# simulate
y = mapply(rbeta,n=1,shape1=betashape[,1],shape2=betashape[,2])
if(keepx) {
return(data.frame(x=x,y=y))
} else {
return(y)
}
}We define function nesufit, which obtains \(\hat{{\theta}}\) fitting the \((X,Y)\) data to the model proposed. By default, the numerical optimization uses a marginal beta estimate as an initial guess; otherwise, supply argument init with your own guess to help the computation.
# function: estimate parameters
nesufit = function(x,y,init=NULL) {
# check compliance to model
nesucheck(theta=init,x=x,y=y)
# encode linear constraints
theui = rbind(c(1,0,0,0),c(-1,1,0,0),c(0,-1,0,0),c(0,0,1,0),c(0,0,-1,1),
c(0,0,0,-1))
theci = c(0,0,-1,0,0,-1)
# get marginal beta as initial guess
if(is.null(init)) {
marginal = MASS::fitdistr(y,densfun="beta",
start=list(shape1=1,shape2=1))$estimate
init = qbeta(c(0.005,0.995,0.005,0.995),shape1=marginal[1],
shape2=marginal[2])
}
# estimate parameters
ans = constrOptim(init,grad=NULL,f=function(theta) {
betashape = nesugetbeta(theta,x)
lik = dbeta(y,shape1=betashape[,1],shape2=betashape[,2])
-1*sum(log(lik))
},ui=theui,ci=theci)
return(ans)
}We define function nesumakeempdist, which creates an empirical sampling distribution for \(\hat{{\theta}}\). Supply argument numreps as the number of replicates. As an initial guess for each replicate, it is best to supply the argument init the estimate obtained by function nesufit.
# function: boostrap to an empirical distribution of parameter estimates
nesumakeempdist = function(x,y,init,numreps) {
# check compliance to model
nesucheck(theta=init,x=x,y=y)
# resample many times
origdat = data.frame(x=x,y=y)
sim = lapply(1:numreps,function(i) {
samp = sample(1:nrow(origdat),size=nrow(origdat),replace=TRUE)
data.frame(x=x,y=y)[samp,]
})
# estimate parameters many times
est = sapply(1:numreps,function(i) {
thetry = try(nesufit(sim[[i]]$x,sim[[i]]$y,init=init)$par)
if(class(thetry)=="try-error") {
message("could not fit at replicate ",i,"; returning ",
"resample instead")
return(sim[[i]])
} else {
thetry
}
})
# collect estimates
ans = data.frame(a0=est[1,],b0=est[2,],a1=est[3,],b1=est[4,])
return(ans)
}We define the function nesutestminhyp, which based on an empirical distribution of \(\hat{{\theta}}\), reports test quantiles relevant to the minimal hypothesis. It is best to supply the empirical distribution obtained by function nesumakeempdist.
# function: verdict on minimal hypothesis
nesutestminhyp = function(a0,b0,a1,b1) {
# check compliance to model
apply(cbind(a0,b0,a1,b1),1,function(theta) {nesucheck(theta,NULL,NULL)})
# get quantile ranks
type1risk = 0.05
quantilerank = type1risk/2 # bonferroni correction
# get test quantiles
floorslope = a1-a0
slopedifference = (b1-b0)-(a1-a0)
test1 = quantile(floorslope,quantilerank)
test2 = quantile(slopedifference,quantilerank)
# report
message(names(test1)," quantile for floor slope is ",test1)
message(names(test2)," quantile for slope difference is ",test2)
return(invisible(NULL))
}We define the function nesuplotpi, which, based on \({\theta}\), draws prediction interval limits over an existing \((U,Y)\) scatterplot. Supply argument yran with the natural range of the data.
# function: draw prediction intervals
nesuplotpi = function(theta,yran=0:1,...) {
# check compliance to model
nesucheck(theta=theta,x=NULL,y=NULL)
theta = setNames(as.list(theta),c("a0","b0","a1","b1")) # assign names
# draw lines
with(theta,lines(0:1,c(b0,b1)*(yran[2]-yran[1])+yran[1],...))
with(theta,lines(0:1,c(a0,a1)*(yran[2]-yran[1])+yran[1],...))
# no return
return(invisible(NULL))
}We define the function nesuplotmean, which, based on \({\theta}\), draws the running mean over an existing \((U,Y)\) scatterplot. Supply argument yran with the natural range of the data.
# function: draw mean
nesuplotmean = function(theta,yran=0:1,...) {
# check compliance to model
nesucheck(theta=theta,x=NULL,y=NULL)
theta = setNames(as.list(theta),c("a0","b0","a1","b1")) # assign names
# function to get mean, for u
getmean = function(u) {
au = with(theta,a0+(a1-a0)*u)
bu = with(theta,b0+(b1-b0)*u)
betashape = nesurb(au,bu)
betashape[,1]/(betashape[,1]+betashape[,2])
}
# draw curve
curve(getmean(x)*(yran[2]-yran[1])+yran[1],from=0,to=1,add=TRUE,...)
# no return
return(invisible(NULL))
}We define the function nesuoverlay, which consolidates the actions for plotting based on the obtained \(\hat{{\theta}}\). Supply argument yran the natural range of the creativity variable. In the scatterplot produced, we follow the legend, to wit: red circles denote actual observations; blue crosses denote simulated observations; solid lines denote prediction interval limits; and dashed lines or curves denote the mean. Note that this function is only for plotting—the simulated data must come from an output of function nesusimulate.
# function: overlay plotting
nesuoverlay = function(theta,actx,acty,simx,simy,yran=0:1,ylab="y",xlab="x") {
# check compliance to model
nesucheck(theta=theta,x=actx,y=acty)
nesucheck(theta=NULL,x=simx,y=simy)
# transform to natural range
untransy = function(y) {
y*(yran[2]-yran[1])+yran[1]
}
# transform to quantile rank
getu = function(x) {
pnorm(x,mean=100,sd=15)
}
# plot background
thexlim = range(c(actx,simx))
theylim = untransy(range(c(acty,simy)))
plot(NA,NA,type="n",xlim=0:1,ylim=theylim,ylab=ylab,xlab=xlab)
# plot points
points(getu(simx),untransy(simy),col="blue",pch=4)
points(getu(actx),untransy(acty),col="red",pch=1)
# plot lines and curves
nesuplotpi(theta=theta,yran=yran,lty=1,lwd=2)
nesuplotpi(theta=theta,yran=yran,lty=1,lwd=2)
nesuplotmean(theta=theta,yran=yran,lty=2,lwd=2)
# no return
return(invisible(NULL))
}The official R package for Necessary Condition Analysis makes it difficult to control graphical settings and output in an interactive session. Instead, we implement the required functionality in our own function nca.
nca = function(x,y,cefdh=TRUE,crfdh=TRUE,draw=FALSE,...) {
# assert
stopifnot(length(x)==length(y))
stopifnot(cefdh|crfdh)
# make cumulative maximum table
sux = sort(unique(x))
cmtab = data.frame(x=sux,y=sapply(sux,function(xval) {
max(y[x<=xval])
}))
cmtab = subset(cmtab,c(TRUE,diff(y)>0))
if(length(cmtab$y)==1) {
message("ceiling zone has zero area; no lines drawn")
return(list(cefdh=0,crfdh=0))
} # for zero-area ceiling scenario
# get scope edges
ymax = max(y)
ymin = min(y)
xmax = max(x)
xmin = min(x)
# get scope area
scopearea = (xmax-xmin)*(ymax-ymin)
# draw scope borders
if(draw) {
lines(c(xmin,xmin,xmax,xmax,xmin),c(ymin,ymax,ymax,ymin,ymin),
lty=3)
}
# CE-FDH
if(cefdh) {
# loop to connect the dots
eceilarea = 0
for(i in 2:nrow(cmtab)) {
eceilarea = eceilarea+(max(y)-cmtab$y[i-1])*
(cmtab$x[i]-cmtab$x[i-1])
if(draw) {
with(cmtab,lines(c(x[i-1],x[i],x[i]),
c(y[i-1],y[i-1],y[i]),...))
}
}
eeffsize = eceilarea/scopearea
}
# CR-FDH
if(crfdh) {
# characterize ceiling
rceilmodcoef = lm(y~x,data=cmtab)$coef
rceil = function(x) {
rceilmodcoef[1]+rceilmodcoef[2]*x
}
rceilinv = function(y) {
(y-rceilmodcoef[1])/rceilmodcoef[2]
}
# identify dots to connect
rconnect = with(new.env(),{
newx = c(xmin,xmax,rceilinv(ymin),rceilinv(ymax))
newx = ifelse(newx<xmin,xmin,ifelse(newx<xmax,newx,xmax))
newy = c(rceil(xmin),rceil(xmax),ymin,ymax)
newy = ifelse(newy<ymin,ymin,ifelse(newy<ymax,newy,ymax))
newdf = data.frame(x=newx,y=newy)
newdf = newdf[with(newdf,order(x,y)),]
})
# loop to connect the dots
rfullarea = 0
for(i in 2:4) {
rfullarea = with(rconnect,rfullarea+
(x[i]-x[i-1])*((y[i]+y[i-1])/2-ymin))
if(draw) {
with(rconnect,{
if(x[i]!=x[i-1]&y[i]-y[i-1]) {
lines(c(x[i-1],x[i]),
c(y[i-1],y[i]),...)
}
})
}
}
reffsize = 1-rfullarea/scopearea
}
# put together
ans = list()
if(cefdh) {
ans = c(ans,list(cefdh=eeffsize))
}
if(crfdh) {
ans = c(ans,list(crfdh=reffsize))
}
return(ans)
}Actual work
With the preliminaries in place, the work in the paper can be reproduced. The first stage is the example scenarios under the model proposed. The second stage is the reanalysis of Jauk et al 2013.
Demonstrative simulations
We produce an ideal triangular scatterplot, \({\theta}= (0.005,0.01,0.01,0.995)\). Result from NCA and segmented regression included.
# ideal triangular
set.seed(2018021101)
datideal = nesusimulate(theta=c(0.005,0.01,0.01,0.995),x=rnorm(500,mean=100,sd=15),keepx=TRUE)
with(datideal,{
plot(x,y,col="blue")
nca(x,y,draw=TRUE)
})## $cefdh
## [1] 0.5366724
##
## $crfdh
## [1] 0.5313165fitseg = with(datideal,segmented::segmented(glm(y~x),seg.Z=~x,psi=median(x)))
fitseg## Call: segmented.glm(obj = glm(y ~ x), seg.Z = ~x, psi = median(x))
##
## Meaningful coefficients of the linear terms:
## (Intercept) x U1.x
## -0.156923 0.002466 0.007098
##
## Estimated Break-Point(s):
## psi1.x
## 90.88
##
## Degrees of Freedom: 499 Total (i.e. Null); 496 Residual
## Null Deviance: 15.29
## Residual Deviance: 8.61 AIC: -601.9plot(fitseg,add=TRUE,rug=FALSE,lwd=2,lty=2,col="black")
We produce the first counterexample, the positive linear relationship, \({\theta}= (0.005,0.3,0.7,0.995)\). Result from NCA included.
# positive linear
set.seed(2018021102)
datcounter1 = nesusimulate(theta=c(0.005,0.3,0.7,0.995),x=rnorm(500,mean=100,sd=15),keepx=TRUE)
with(datcounter1,{
plot(x,y,col="blue")
nca(x,y,draw=TRUE)
})
## $cefdh
## [1] 0.4064131
##
## $crfdh
## [1] 0.4034443We produce the second counterexample, the right-opening megaphone, \({\theta}= (0.4,0.6,0.005,0.995)\). Result from NCA included.
# megaphone
set.seed(2018021103)
datcounter2 = nesusimulate(theta=c(0.4,0.6,0.005,0.995),x=rnorm(500,mean=100,sd=15),keepx=TRUE)
with(datcounter2,{
plot(x,y,col="blue")
nca(x,y,draw=TRUE)
})
## $cefdh
## [1] 0.2168819
##
## $crfdh
## [1] 0.204522We produce the third counterexample, the less-than-perfect triangular shape, \({\theta}= (0.005,0.1,0.15,0.995)\). Result from NCA included.
# unideal triangular
set.seed(2018021104)
datcounter3 = nesusimulate(theta=c(0.005,0.1,0.15,0.995),x=rnorm(500,mean=100,sd=15),keepx=TRUE)
with(datcounter3,{
plot(x,y,col="blue")
nca(x,y,draw=TRUE)
})
## $cefdh
## [1] 0.4262973
##
## $crfdh
## [1] 0.4210306Reanalysis of Jauk et al 2013
We load the data from Jauk et al 2013. We exclude observations that have missing IQ scores.
# load data
dato = read.csv("./fromjauketal2013.csv",header=TRUE) # raw
dat = subset(dato,!is.na(IQ)) # complete observations onlyOf the data, we take descriptives. We also show the scatterplots and results with NCA.
# get descriptives
t(sapply(dat,function(u) {
list(mean=mean(u),sd=sd(u),coriq=cor(u,dat$IQ),min=min(u),max=max(u))
}))## mean sd coriq min max
## IQ 107.2142 14.62614 1 59.27333 147.3661
## CP_avgOriginality 1.823315 0.1214644 0.3521672 1.451735 2.162293
## CA 40.72054 35.15121 0.2814005 0 208# raw scatterplots and NCA
with(dat, {
plot(IQ,CP_avgOriginality,xlab="IQ",ylab="Creative Potential")
nca(IQ,CP_avgOriginality,draw=TRUE)
})
## $cefdh
## [1] 0.1880758
##
## $crfdh
## [1] 0.1735106with(dat, {
plot(IQ,CA,xlab="IQ",ylab="Creative Achievement")
nca(IQ,CA,draw=TRUE)
})
## $cefdh
## [1] 0.3889667
##
## $crfdh
## [1] 0.3380747The same procedure was done separately for IQ vs. CPA and IQ vs. CA. The general workflow is as follows.
- The creativity variable is transformed in line with the model proposed, i.e. \(0 \leq Y \leq 1\).
- Given \((X,Y)\) data, \(\hat{{\theta}}\) is obtained.
- To visually inspect the model fit, \({Y \rvert \{X = x\}}\) was simulated from the value of \(\hat{{\theta}}\) and a resampling of actual data on \(X\) plus uniform jitter. Two scatterplots are produced, one in terms of raw IQ \(X\) and the other in terms of IQ quantile rank \(U\). Each scatterplot shows actual data, simulated data, prediction intervals, and the conditional expectation.
- By bootstrap, an empirical distribution of \(\hat{{\theta}}\) was obtained. This empirical distribution is tested with respect to the minimal hypothesis.
For IQ vs. CPA, the minimal hypothesis is disconfirmed.
# IQ vs. CPA
Sys.time() # time start## [1] "2018-02-22 11:50:44 +08"datcpa = with(dat,{data.frame(x=IQ,
y=(CP_avgOriginality-1)/(4-1))}) # set up and transform
resucpa1 = with(datcpa,nesufit(x,y)) # estimate## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs producedresucpa1## $par
## [1] 0.1416158 0.3620695 0.2072610 0.3894150
##
## $value
## [1] -551.6706
##
## $counts
## function gradient
## 229 NA
##
## $convergence
## [1] 0
##
## $message
## NULL
##
## $outer.iterations
## [1] 2
##
## $barrier.value
## [1] 0.0001834512set.seed(2018021111)
resucpa2 = nesusimulate(theta=resucpa1$par,
x=sample(datcpa$x,size=500,replace=TRUE)+runif(500,-2,+2),keepx=TRUE) # simulate
nesuoverlay(theta=resucpa1$par,actx=datcpa$x,acty=datcpa$y,simx=resucpa2$x,
simy=resucpa2$y,yran=c(1,4),ylab="Creative Potential",
xlab="IQ") # plot simulation
set.seed(2018021112)
resucpa3 = nesumakeempdist(x=datcpa$x,y=datcpa$y,numreps=400,
init=resucpa1$par) # bootstrap
with(resucpa3,nesutestminhyp(a0,b0,a1,b1)) # verdict## 2.5% quantile for floor slope is 0.0383664570873507## 2.5% quantile for slope difference is -0.0753592139472031Sys.time() # time end## [1] "2018-02-22 16:17:03 +08"For IQ vs. CA, the minimal hypothesis is confirmed. Several observations were excluded due to having a CA score of zero. Note the time it takes to produce an empirical distribution of 400 replicates.
# IQ vs. CA
Sys.time() # time start## [1] "2018-02-22 16:17:05 +08"datca = with(subset(dat,CA>0),{data.frame(x=IQ,
y=(CA-0)/(440-0))}) # set up and transform
resuca1 = with(datca,nesufit(x,y))## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs producedresuca1 # estimate## $par
## [1] 0.000476356 0.207292768 0.001956318 0.455309778
##
## $value
## [1] -421.4422
##
## $counts
## function gradient
## 340 NA
##
## $convergence
## [1] 0
##
## $message
## NULL
##
## $outer.iterations
## [1] 3
##
## $barrier.value
## [1] 0.0001215479set.seed(2018021121)
resuca2 = nesusimulate(theta=resuca1$par,
x=sample(datca$x,size=500,replace=TRUE)+runif(500,-2,+2),keepx=TRUE) # simulate
nesuoverlay(theta=resuca1$par,actx=datca$x,acty=datca$y,simx=resuca2$x,
simy=resuca2$y,yran=c(0,440),ylab="Creative Achievement",
xlab="IQ") # plot simulation
set.seed(2018021122)
resuca3 = nesumakeempdist(x=datca$x,y=datca$y,numreps=400,
init=resuca1$par) # bootstrap
with(resuca3,nesutestminhyp(a0,b0,a1,b1)) # verdict## 2.5% quantile for floor slope is 0.00110834226870441## 2.5% quantile for slope difference is 0.132477313217582Sys.time() # time end## [1] "2018-02-22 17:53:00 +08"End of session
At this point, reproducing the work is complete. We save our work and show session info.
# session end
save.image("./workspace.RData") # save workspace
devtools::session_info() # session info## Session info -------------------------------------------------------------## setting value
## version R version 3.2.3 (2015-12-10)
## system x86_64, linux-gnu
## ui X11
## language en_PH:en
## collate en_PH.UTF-8
## tz Asia/Manila
## date 2018-02-22## Packages -----------------------------------------------------------------## package * version date source
## backports 1.1.2 2017-12-13 CRAN (R 3.2.3)
## base * 3.2.3 2016-01-14 local
## datasets * 3.2.3 2016-01-14 local
## devtools * 1.13.4 2017-11-09 CRAN (R 3.2.3)
## digest 0.6.14 2018-01-14 CRAN (R 3.2.3)
## evaluate 0.10.1 2017-06-24 CRAN (R 3.2.3)
## graphics * 3.2.3 2016-01-14 local
## grDevices * 3.2.3 2016-01-14 local
## htmltools 0.3.6 2017-04-28 CRAN (R 3.2.3)
## knitr 1.18 2017-12-27 CRAN (R 3.2.3)
## magrittr 1.5 2014-11-22 CRAN (R 3.2.3)
## MASS * 7.3-45 2015-11-10 CRAN (R 3.2.2)
## memoise 1.1.0 2017-04-21 CRAN (R 3.2.3)
## methods * 3.2.3 2016-01-14 local
## Rcpp * 0.12.15 2018-01-20 CRAN (R 3.2.3)
## rmarkdown 1.8 2017-11-17 CRAN (R 3.2.3)
## rprojroot 1.3-2 2018-01-03 CRAN (R 3.2.3)
## segmented * 0.5-3.0 2017-11-30 CRAN (R 3.2.3)
## stats * 3.2.3 2016-01-14 local
## stringi 1.1.6 2017-11-17 CRAN (R 3.2.3)
## stringr 1.2.0 2017-02-18 CRAN (R 3.2.3)
## tools 3.2.3 2016-01-14 local
## utils * 3.2.3 2016-01-14 local
## withr 2.1.1 2017-12-19 CRAN (R 3.2.3)
## yaml 2.1.16 2017-12-12 CRAN (R 3.2.3)