#_____________________________________________________________________________________________
# R code accompanying Verkuilen, Smithson, Kievit, Zand Scholten
#_____________________________________________________________________________________________
# This simulation illustrates the assignment of membership values to a fuzzy set in the
# context of diagnosis of Major Depressive Disorder (MDD).
# First, np persons with different probabilities of displaying ns symptoms are simulated.
# Second, these probabilities are transformed to 'observed' scores indicating whether a
# symptom is actually displayed or not.
# Third, membership to the set of clinically depressed persons (displaying MDD) is assigned.
# Assigment is performed using different methods:
# 1) DSM procedure: person is classified as depressed when displaying 5 or more symptoms
# 2) Direct Assignment: person is assigned value by summing weighted symptom scores
# 3) Assignment using indirect scaling of symptoms (iss): This is the weighting summation of the 9 symptoms.
#_____________________________________________________________________________________________
# SIMULATION OF PERSONS, SYMPTOMS, SYMPTOM PROBABILITIES AND SYMPTOM SCORES
np = 10000; # Number of simulated persons
ns = 9; # Number of symptoms
a = rep(1,ns); # Discrimination parameters used to generate symptom probabilities
b = c(-1.80,-0.58,-0.37,
-0.21,-0.21,-0.18,
0.90, 0.95,2.47);# Difficulty parameters used to generate symptom probabilities (see Aggen et al.)
t = rchisq(np,2); # Chisquare distributed latent values
t = ((6*t)/max(t)-3); # Rescaling to range [-3,3]
prob = matrix(,np,ns); # Matrix will contain probabilities of displaying each of ns symptoms for np persons
# Generates probability of displaying a certain symptom according to 2pl IRT model for each of np persons
for(s in 1:ns) { # For each symptom s
prob[,s]=exp(a[s]*(t-b[s]))/(1+exp(a[s]*(t-b[s])));
}
score = matrix(,np,ns) # Matrix will contain dichotomized scores, indicating whether symptom is displayed (1) or not (0)
# Generates dichotomous score by sampling from a binomial distribution for each of np persons with probability determined earlier
for(s in 1:ns) { # For each symptom s
score[,s]=rbinom(np,1,prob[,s]);
}
# Symptom names
colnames(score)=c("Depressed", "Weightapp", "Nointerest", "Sleepprob", "Fatigue", "Psycmotor", "Worthless", "Concentra", "Suicidal");
# All data aggregated in a data frame
data = data.frame(cbind('Theta'=t, score));
#_____________________________________________________________________________________________
# DIFFERENT ASSIGMENT METHODS
#1 Assigment using DSM criteria (dsm). Assuming these symptoms are all longer than two weeks, >5=depressed
sms = rowSums(score); # summing into a symptom count for each person
dsm = ifelse(sms>5,1,0); # recoding into 1 if symptom count > 5, else 0
#2 Direct assignment (das): Symptoms are weighted by using linear regression weights obtained from a sample of persons who
# are assigned membership values directly by experts. Here we use the original latent scores as direct assignment values
data = data.frame(cbind('Theta'=t, score, sms, dsm, iss)); # All data aggregated in a data frame
# beta weights obtained from linear regression of theta on ns symptom scores
bts = lm(Theta ~ Depressed + Weightapp + Nointerest + Sleepprob + Fatigue + Psycmotor + Worthless + Concentra + Suicidal, data)$co;
bwts = matrix(rep(bts[-1],np),np,ns,T); # Matrix of beta weights repeated np times
bwsc = bwts*score; # Symptom scores for np persons multiplied by beta-weights
das = rowSums(bwsc); # The weighted sum score
das = das/sum(bts[-1]); # Dividing weighted sumscore by max score, rescaling to range[0,1]
data = data.frame(cbind(data,das)); # Adding das score to data
layout(matrix(1:2,1,2));
hist(t)
hist(rowSums(score))
#3 Assignment using indirect scaling of symptoms (iss): This is the weighting summation of the 9 symptoms.
#weights = matrix(rep(c(.03,.04,.07,.09,.11,.12,.13,.17,.24),np),np,ns,T); # Symptoms weights that sum to 1. These are example weights and may be adapted
weights = matrix(rep((b+1.8)/sum(b+1.8),np),np,ns,T); # Weights based on Aggen et al.
wscore = weights*score; # multiplying weights by presence of symptom for each person to get weigted score
iss = rowSums(wscore); # the weighted sum score