#_____________________________________________________________________________________________
# 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