利用 GAMLSS 对心理测验进行基于回归的常模分析

文摘   2024-10-16 08:48   中国  


PSYCH统计实验室


前言

标准化的心理测试被广泛使用。它们在个体评估和心理研究中起着至关重要的作用。高质量标准化测试的一个核心特点是它们具有良好的常模分数。大多数心理测试都有常模参照分数,包括智力测试、发展测试、神经心理测试以及临床测试。在常模参照测试中,原始分数被转换为常模参照分数,这表达了个体的表现相对于测试的参照人群的表现。

常模参照心理测试的一个显著特点是,这些常模通常依赖于年龄(例如,智力测试)有时还依赖于性别和/或教育水平(例如,神经心理学测试)。这意味着实际上存在多个需要常模化的参照人群,它们共同构成了测试的常模人群。这些常模是在测试构建阶段建立的,基于在常模人群样本中收集的分数。

在常模参照测试中的常模可能依赖于个体特征,如年龄和性别,研究将这些常模所依赖的变量称为常模预测变量。近年来,连续常模已被测试构建者采纳,用于计算他们的常模,因为它在准确性方面具有有利的特性。连续常模的关键是,在计算常模时,明确使用了常模预测变量(如年龄或教育水平)的连续(或有序分类)性质所提供的信息。即使涉及一个或多个离散常模预测变量,也可以使用连续常模。即使对于单一离散常模预测变量(例如,仅性别),也可以应用连续常模,但与传统常模相比,其优势会变得小得多。

基于回归的常模分析相较于以往的方法(如半参数常模分析和推断常模分析)有两个显著优势:首先,它有现成的统计标准用于模型选择和模型评估。其次,在计算过程中,没有必要将连续变量如年龄离散化,从而避免了对区间宽度的任意且可能有影响的决定。而利用GAMLSS (广义可加模型位置、尺度和形状)进行基于回归的常模分析则几乎适用于广泛的实证常模分析案例,并包含了迄今为止在基于回归的常模分析中使用的所有特定模型,因为GAMLSS可以实现针对不同的分布模式(如BB分布、BCPE分布等)的线性/非线性建模。

那么如何利用GAMLSS进行基于回归的常模分析呢?Timmerman等人总结出了六个标准步骤,如下表所示。


步骤详解


1

第一步:定义测验的参考人群、

常模人群和目标人群

常模分数是相对于测试的参考人群来表达的。例如,参考人群可以是与受测者同龄的国家一般人群。所有参考人群共同(例如,测试目标人群年龄范围内的所有年龄)构成了测试的常模人群。测试的常模人群(即,测试的所有参考人群共同)可能等于目标人群,或者可能是目标人群的一部分。例如,针对社区人群的智力测试具有相同的目标人群和常模人群。神经心理学测试通常以健康个体作为常模人群,而测试的目标人群包括健康和不健康的个体。

2

第二步:设计研究以收集常模

样本数据

一个良好构成的常模样本对于获得正确的常模至关重要。正确的常模是无偏的且足够精确的。这些只有在无偏且足够精确估计的原始测试分数的分布条件下才能实现,这些分布条件是基于常模预测变量的。众所周知,从总体中随机抽样可以实现无偏估计,但现实是随机抽样不可能实现,因此,某些情况下,分层抽样是一个良好的替代方案,但由于非相应偏差的存在,也会造成偏差。一种实际上可行的随机抽样替代方法是判断抽样

3

第三步:选择一个候选的

GAMLSS分布作为条件原始

分数模型 

基于回归的常模分析的核心步骤是将原始测试分数分布建模为常模预测变量函数的适当模型。在常模分析中所需的建模类型与标准(线性)回归建模类似,但有两个重要区别。一个是不一定需要服从正态分布,另一个是既可以建模线性关系,也可以建模非线性关系,因此,这需要使用比标准线性回归模型更灵活的模型。

为了在常模样本中适当地模拟原始测试分数分布,需要一种灵活的建模方法,这在GAMLSS中可以找到。GAMLSS是单变量分布回归模型,其中假定分布的所有参数都可以作为预测变量的加性函数来建模。例如,当使用正态分布时,可以将均值(μ)和标准差(σ)都作为预测变量的函数来建模。GAMLSS可以在R(中使用gamlss包进行拟合。该包包括超过100种连续的、离散的和混合类型的分布,用于建模结果变量。这提供了很大的建模灵活性,使得找到合适的模型变得相当可能。

要指定一个GAMLSS模型,需要选择一个可能的拟合分布。通过与原始测试分数的性质匹配,考虑分数范围以及它们的离散或连续特性,可以确定适合使用的分布。通常,选择尽可能少参数的分布是明智的,以避免过度拟合。此外,一旦选择了GAMLSS分布,就需要为每个分布参数确定所谓的链接函数。链接函数与分布函数共同决定了参数值的可能范围。非常流行的函数是恒等链接,它保留了参数值在其原始范围内,以及对数链接函数,它将参数值限制为非负。通常,恒等链接函数用于位置相关的参数(如μ)和偏度(如σ),对数链接函数用于尺度和峰度。

一种常见出现的原始量表类型是有序分类量表,具有固定的最小和最大分数。当测试分数是构成量表的项目分数之和时,就会出现这种类型。典型的项目分数是二元的(例如,0 表示错误;1 表示正确),或者是在李克特量表中的多个有序类别。对于受预测变量限制的分数范围(例如,最大25分),贝塔二项分布(BB)对于许多有序分类量表来说,拟合得相当好。如果一个量表由二元项目组成,其中个体能力的分布遵循贝塔分布,项目难度相等,或者构成量表的项目难度遵循贝塔分布,个体能力相等,那么贝塔二项分布是有理论依据的。在这两种假设下,观察到的量表分数遵循贝塔二项分布。在由二元项目组成的量表的标准化模型中,项目难度和个体能力水平可能会在项目和个体之间变化。由于条件化了常模预测变量,能力水平的变异性通常会受到一定程度的限制。此外,即使这些假设没有得到满足,贝塔二项分布也可能提供适当的统计拟合。然而,评估模型拟合仍然很重要。

连续分布适用于模拟在连续尺度上测量的量表,以及近似具有足够多类别的分类尺度(例如,超过25个类别)。Box-Cox幂指数(BCPE)分布是一种灵活的、实践中通常拟合度很好的连续分布。例如,BCPE已被用于模拟涉及反应时间(RT)作为神经心理学测试COTAPP测试分数的量表。使用BCPE分布的理由是统计性质的,因为它是一种灵活的分布,可以拟合广泛的连续经验分布。BCPE分布有四个参数,与位置(μ,中位数)、尺度(σ,近似变异系数)、偏度(γ,对称变换)和峰度(δ,幂指数参数)有关。请注意,在BCPE分布的上下文中,符号μ和σ指的是中位数和尺度,因此它们的含义与默认值(即均值和标准差)不同。正态分布是BCPE分布的一个特殊情况,即当γ = 1且δ = 2时。因此,BCPE分布是连续尺度和具有足够多类别的分类尺度的首选

4

第四步:选择候选函数将常模

预测变量与GAMLSS分布参数

相关联

对于一个连续的预测变量,一个简单的方法是使用线性模型来拟合每个参数。例如,使用正态分布将导致参数μ和σ都线性依赖于预测变量。注意,如果σ被建模为截距(即,视为与任何常模预测无关),模型将简化为标准的线性回归模型。然而常模预测中的非线性关系相当常见,这使得需要对非线性进行建模。

涉及常模预测变量和结果变量之间平滑关系的非线性,可以使用多项式样条函数来建模。多项式是最简单的方法,涉及在线性方程中添加一个或多个高阶项(例如,年龄平方age^2,年龄立方age^3)。为了避免由于预测因子的多重共线性导致的估计问题,建议使用预测因子的中心化版本和/或常模预测因子集的正交化版本。

样条函数是分段多项式函数,用于转换常模预测变量。在回归中使用这些转换后的常模预测变量可以得到一个平滑的估计函数。这样的函数可以采取任何平滑的函数形式。使用样条函数的关键问题在于所需的平滑度,以实现一个能很好地代表总体的模型。因此,与任何模型一样,需要平衡欠拟合和过拟合。不同的样条类型以不同的方式操纵平滑度。一种流行的是P样条,因为它具有有利的属性(例如,数值稳定且易于实现)。此外,它只需要一个惩罚参数来操纵整个函数的平滑度,使得P样条易于应用。如果事先知道平滑函数是单调递增(或递减)的,那么单调P样条变体可以用来实现更有效的估计。当建模已知随年龄逐渐增加的测试分数时,例如儿童时期的智力分数,这种单调性约束特别有用。在GAMLSS模型中,这可以通过将位置参数(μ)作为年龄的函数来表达单调性约束。请注意,对其他参数如散布和偏度施加单调性约束通常是不合适的,因为这些通常不显示单调模式。

5

第五步:执行模型选择以得到

估计的GAMLSS模型

使用多项式或P样条拟合函数都需要模型选择。这归结为选择多项式回归中包含的多项式,或者在P样条中选择惩罚参数的值。目前还没有关于如何选择这些参数的普遍共识。一种流行的方法是从一系列可能拟合良好的模型中,使用统计标准选择一个或两个有前景的候选模型。在这些候选模型中,然后通过视觉诊断来评估模型拟合情况。

存在两种类型的统计模型选择标准。两者都旨在适当地拟合样本数据,同时防止过拟合,以确保模型对总体的泛化能力。在交叉验证中,这是通过直接评估模型对新数据的预测质量来完成的。在广义赤池信息准则(GAIC)中,这是通过间接地惩罚模型中包含的参数数量来完成的。不同的GAIC变体在惩罚程度上有所不同,例如赤池信息准则(AIC)、GAIC(3)和贝叶斯信息准则(BIC),对参数数量的惩罚分别等于2、3和ln(N)。研究者推荐使用贝叶斯信息准则(BIC)作为主要的统计标准来识别候选模型,并使用另一个广义信息准则(例如,AIC)来评估所青睐模型的一致性。

6

第六步:根据估计的GAMLSS

模型计算量表的常模分数

将原始测试分数转换为任何所需的常模分数是基于累积分布函数(CDF)进行的。从估计的GAMLSS模型中,可以为每个预测值(或者在多个预测值的情况下,任何预测值的组合)获得一个模型隐含的CDF。因此,可以为任何所需的参考人群推导出常模。请注意,不鼓励在观察数据的边界之外进行外推。

操作


###所有代码均来源于(Timmerman et al, 2021)研究

#首先导入3个包,如果未安装请先安装

library(gamlss)

library(gamlss.tr)

library(scales)

#定义str_eval的函数,其作用是将字符串形式的表达式作为数学表达式进行求值

str_eval=function(x) {return(eval(parse(text=x)))}

#设置工作目录

setwd("D:\\Rdaima\\Regression-Based Norming\\met0000348_Supplemental-Material")

#导入数据

IDS2_sample <- read.table("IDS2_sample.txt", header = T)

#数据框IDS2_sample如下所示


#将y14与年龄的关系用散点图绘制出来

plot(IDS2_sample$y14, x = IDS2_sample$age, pch = 16, col = scales::alpha("black", 0.6), ylab = "Test score subtest 14", xlab = "Age", ylim = c(0, 34))


#######################拟合BB分布 ########################

max_score <- 34 # 定义理论最大值,y14的理论最大值是34

IDS2_sample$BB_y14 <- cbind(IDS2_sample$y14, max_score-IDS2_sample$y14) #定义适用于#BB分布的结局指标格式

# 比较mu和sigma的所有与年龄无关的模型组合,以及mu和sigma的四次正交多项式。max_mu_degree <- 4

max_sigma_degree <- 4

i <- 0

namesBB <- expand.grid(sigma = 0:max_sigma_degree, mu = 0:max_mu_degree)#排列组合

# 创建空列表存储模型.

mod_BB_list <- sapply(paste(namesBB$mu, namesBB$sigma, sep = ""), function(x) NULL)

for (d_mu in 0:max_mu_degree){

  if (d_mu == 0){

    mod_mu_degree <- "1"

  } else {

    mod_mu_degree <- paste("poly(age, ", d_mu,")", sep = "")

  }

  for (d_sigma in 0:max_sigma_degree){

    if (d_sigma == 0){

      mod_sigma_degree <- "1"

    } else {

      mod_sigma_degree <- paste("poly(age, ", d_sigma,")", sep = "")

    }

    tryCatch({ # START TRYCATCH (If a model cannot be estimated, continue with the next model)

    i <- i + 1

    mod_BB_list[[i]] <- str_eval(paste("gamlss(BB_y14 ~ ",mod_mu_degree,", sigma.formula = ~ ",mod_sigma_degree,", family = 'BB', data = IDS2_sample, method = RS(1000), pb = max_score)", sep = ""))

    }, warning= function(w) Inf) ### END TRYCATCH

    }

}

mod_BB_list <- mod_BB_list[!sapply(mod_BB_list, is.null)] # 删除所有NULL元素(无法估计的模型)

BIC_BB <- unlist(lapply(mod_BB_list, BIC)) # 从每个模型中提取BIC

BB_mod <- mod_BB_list[[which(BIC_BB == min(BIC_BB))]] # 选择最低BIC的bb模型

 

# 还为mu和sigma拟合了一个具有P样条的模型;μ的P样条单调递增:“pbm(…,mono=”up“)”,使用BIC选择惩罚

BB_mod_sp <- gamlss(BB_y14 ~ pbm(age, method = "GAIC", k = log(nrow(IDS2_sample)), inter = 5, mono = "up"), sigma.formula = ~ pb(age, method = "GAIC", k = log(nrow(IDS2_sample)), inter = 5), family = "BB", data = IDS2_sample, method = RS(1000), pb = max_score)

 

################## 拟合BCPE分布##################

 

IDS2_sample$y14_a <- IDS2_sample$y14 + 0.0001 # BCPE无法处理等于0的分数。为了估计测试分数0的百分位数,我们创建了一个调整后的结果变量(“y14_a”),其中所有测试分数都增加了0.0001。这样,模型中原始测试分数0.0001的估计百分位数等于原始测试分数0的估计百分数。

 

### Applied free order model selection procedure to y14_a -> see separate R file【这一块使用free order model selection procedure进行最佳的BCPE模型参数选择,因为如果继续使用for循环效率太低,这块相关代码可参考Timmerman等人的研究,这里不给出了】

 

# 所选模型:2-1-0-0(=μ的二阶年龄多项式,σ的一阶年龄多项式(=线性效应),nu和tau的截距)

BCPE_mod <- gamlss(y14_a ~ poly(age, 2),

                   sigma.formula = ~ poly(age, 1),

                   nu.formula =~ 1,

                   tau.formula =~ 1,

                   family = "BCPE", data = IDS2_sample, method = RS(1000))

 

# 用μ的P样条而不是二阶多项式来估计BCPE分布

BCPE_mod_sp <- gamlss(y14_a ~ pbm(age, method = "GAIC", k = log(nrow(IDS2_sample)), inter = 5, mono = "up"),

                        sigma.formula = ~ age, # = linear effect, like poly(age, 1)

                        nu.formula = ~ 1,

                        tau.formula = ~ 1,

                        family = BCPE, data = IDS2_sample, method = RS(1000))

 

############# 拟合向右截断的BCPE分布 #############

 

BCPEtr34 <- gen.trun(c(0, 34), family = "BCPE", type = "both", name = "tr34") # 生成截断为0和34的BCPE分布。这与34处的右截断完全相同,因为此BCPE分布已经在0处左截断。

 

## Free order model selection with BCPEtr34 distribution -> 2 - 1 - 0 - 0 (same as in "BCPE_mod")

BCPE_mod_tr <- gamlss(y14_a ~ poly(age, 2),

                   sigma.formula = ~ poly(age, 1),

                   nu.formula =~ 1,

                   tau.formula =~ 1,

                   family = BCPEtr34, data = IDS2_sample, method = RS(1000))

 

# 用P样条拟合截断的BCPE分布

BCPE_mod_tr_sp <- gamlss(y14_a ~ pbm(age, method = "GAIC", k = log(nrow(IDS2_sample)), inter = 5, mono = "up"), sigma.formula = ~ age, nu.formula =~ 1, tau.formula =~ 1,family = BCPEtr34, data = IDS2_sample, method = RS(1000))

####################### 拟合标准正态分布#######################

 

# Normal distribution (NO) for normality

NO_mod <- gamlss(y14 ~ age, # linearity

                    sigma.formula = ~ 1, # homoscedasticity

                    family = NO, data = IDS2_sample, method = RS(1000))

 

#########################比较估计的模型#########################

 

### 相对拟合指标

AIC(BB_mod, BB_mod_sp, BCPE_mod, BCPE_mod_sp, BCPE_mod_tr, BCPE_mod_tr_sp, NO_mod) # 比较AIC

BIC(BB_mod, BB_mod_sp, BCPE_mod, BCPE_mod_sp, BCPE_mod_tr, BCPE_mod_tr_sp, NO_mod)# 比较BIC

 

### 绝对拟合指标

wp(BB_mod_sp, xvar = IDS2_sample$age, ylim.worm = 1.5)

title(main='BB_mod_sp', line = -5.2)

#在图中,展示了默认BCPE模型使用P样条(面板B)的蠕虫图。蠕虫图是一系列去趋势的Q-Q图,如果常模预测值(此处:年龄)被分为子范围。蠕虫图可视化了统计模型对数据的拟合程度,用于找出可以改进拟合的位置,以及比较不同模型的拟合情况。蠕虫图上方的蓝色条形表示每个面板中的预测范围,这些范围按行从左下到右上排序。BCPE模型的蠕虫图(面板B)显示大多数观察值(即点)在95%的带状区域内,表明该模型对所有年龄范围的拟合相当好,除了8到12岁的年龄范围(右下角的图)。局部拟合不佳可能是由于在9岁时,观察到的原始测试分数低于10的极端少量。

#####第6步,计算常模分数

#预测5-21范围内1000个预测值的分布参数

pop_age <- data.frame(age = seq(5, 21, length.out = 1000))

pred_distr <- predictAll(BCPE_mod_sp, newdata = pop_age, type = "response")

 

# 估计原始测试分数在0-34范围内的百分位数,步长为1。

# 手动将分数增加0.0001(因为分数0是不可能的)

#因此,百分位数也被评估为分数0.0001-34.0001。

Min_score <- 0.0001 # Minimum (adjusted) raw test score for which to calculate percentiles

Max_score <- 34.0001 # Maximum (adjusted) raw test score for which to calculate percentiles

step_size <- 1

 

#创建空矩阵存储百分数; columns = raw test scores & rows = age values.

CDF_matrix <- matrix(NA, ncol = (Max_score - Min_score)/step_size + 1, nrow = nrow(pop_age))

# 基于年龄的估计分布参数的估计百分位数;

# 使用BCPE分布的分布函数: "pBCPE()".

for (i in 1:nrow(pop_age)){

  CDF_matrix[i,] <- pBCPE(seq(from = Min_score, to = Max_score, by = step_size), mu = pred_distr$mu[i], sigma = pred_distr$sigma[i], nu = pred_distr$nu[i], tau = pred_distr$tau[i], lower.tail = TRUE) # lower.tail = FALSE when a higher test score must correspond to a lower rather than higher percentile.

}

plot(100*CDF_matrix[,1], x = pop_age[,"age"], ylim = c(0, 100), ylab = "Percentile", xlab = "Age", type = "l")

for (i in 2:ncol(CDF_matrix)){

  lines(100*CDF_matrix[,i], x = pop_age[,"age"])

}

### 将样本观测值添加到图中

# 计算标准样本中观测值的百分位数

CDF_sample <- rep(NA, times = nrow(IDS2_sample))

for (i in 1:nrow(IDS2_sample)){

  CDF_sample[i] <- pBCPE(IDS2_sample$y14_a[i], mu = BCPE_mod_sp$mu.fv[i], sigma = BCPE_mod_sp$sigma.fv[i], nu = BCPE_mod_sp$nu.fv[i], tau = BCPE_mod_sp$tau.fv[i], lower.tail = TRUE)

}

# 加观测值到图中

points(100*CDF_sample, x = IDS2_sample$age, col = "aquamarine4", pch = 1)

# 将估计百分位数转换为估计的归一化z分数

z_matrix14 <- qnorm(CDF_matrix)

 

# 将标准化的z分数限制在一定范围内[-3; +3]

z_matrix14 <- apply(z_matrix14, 1, function(x) {ifelse(x > 3, 3, x)})

z_matrix14 <- apply(z_matrix14, 1, function(x) {ifelse(x < -3, -3, x)})

 

# 将样本中的百分位数转换为归一化的z分数; 并且限制到 [-3; +3]

z_sample14 <- ifelse(qnorm(CDF_sample) > 3, yes = 3,

                   no = ifelse(qnorm(CDF_sample) < -3, yes = -3, no = qnorm(CDF_sample)))

z_sample14<-as.data.frame(z_sample14)

# Centile plots

centiles.fan(BCPE_mod_sp, xvar = IDS2_sample$age, points = T, colors = "gray", col = "black", pch = 16, ylab = "Test score", xlab = "Age", main = "Centile curves using BCPE with P-splines")

legend("bottomright", inset = c(0,0), title = "Percentiles", c("25-50th & 50-75th", "10-25th & 75-90th", "2-10th & 90-98th", "0.4-2th & 98-99.6th"),

       pch = c(15, 15), pt.cex = 1.5, col = c("gray30", "gray55", "gray75", "gray95"), cex = 0.7, title.adj = 0.1)

到此为止,我们实现了构建每个特定年龄个体y14的常模分数了。



参考内容

Timmerman, M. E., Voncken, L., & Albers, C. J. (2021). A tutorial on regression-based norming of psychological tests with GAMLSS. Psychological methods, 26(3), 357–373. https://doi.org/10.1037/met0000348


PSYCH统计实验室

通知公告

网络分析课程目前开放视频课啦

单次课200元/讲(学生),250元/讲(非学生)

共有四讲内容:

①横断面网络分析简介与基础

②网络分析与因子分析

③交叉滞后网络分析

④时间序列网络分析

购买后开放视频权限14天,可多次申请。

并赠送所有课程相关资料(无PPT)

如果想申请购买,请联系M18812507626


更多资讯

关注我们

文稿:摘星

排版:Little Star

责编:Wink
审核:摘星

本文由“Psych统计自习室”课题组原创,欢迎转发至朋友圈。如需转载请联系后台,征得作者同意后方可转载。


Psych统计自习室
大家好,我们是由来自北京师范大学,西南大学,天津医科大学等高校在读硕士、博士研究生组成的一个科研团队——Psych统计自习室。Psych统计自习室旨在关注心理学、精神病学领域的最前沿的系列研究,并做前沿统计知识的分享。
 最新文章