# Inspecting mechanism for missing values in categorical data without prior knowledge

## Scenario

I am inspecting the Soybean data set, which has a quite a number of missing values for various categorical variables.

## Plan

My plan is to eventually perform data imputation. However, currently I am trying to understand the mechanism behind these missing values (MAR, MNAR) using graphs, as will be presented below.

## Problem

I feel I have come to a standstill with regards to my exploration: I am unable to decide whether some data is MAR or MNAR after attempting several graphs. I would very much appreciate any advice or insights for understanding the mechanism behind my missing data values.

### Code and graphs

#load required libraries and data set
require(mlbench)
require(ggplot2)
require(vcd)
require(reshape2)

data(Soybean)
d1<- Soybean

#generate barplots for number of missing values and their proportion against variables
var.missing<- sapply(d1,function(x)sum(is.na(x)))
var.missing<- var.missing[order(var.missing)]
missing.df<- data.frame(variable=names(var.missing),missing=var.missing,missing.prop=var.missing/dim(d1)[1],stringsAsFactors=FALSE)
missing.df$variable<- factor(missing.df$variable,levels=missing.df$variable,ordered=FALSE) g1<- ggplot(data=missing.df,aes(x=variable,y=missing)) + geom_bar() + labs(x="Variables",y="Number of missing values") + theme(axis.text.x=element_text(angle=45, hjust=1))  g2<- ggplot(data=missing.df,aes(x=variable,y=missing.prop)) + geom_bar() + labs(x="Variables",y="Proportion of missing values") + theme(axis.text.x=element_text(angle=45, hjust=1)) #so, lodging, seed.tmt, sever, hail have most missing values #then come: germ, leaf.mild, shriveling, seed.discolor, fruiting.bodies, leaf.shred, #seed.size, mold.growth, seed, fruit.pods, lead.malf, leaf.size, leaf.marg, leaf.halo  #let's check the proportion of missing values per class df.per.class<- split(d1,d1$Class)
rows.per.class<- sapply(df.per.class, nrow)
tot.values.per.class<- sapply(rows.per.class,function(x)x*dim(d1)[2])
miss.per.class.usingRows<- sapply(df.per.class,function(x)apply(x,1,function(y)sum(is.na(y))))
miss.rows.per.class<- sapply(miss.per.class.usingRows,function(x)sum(x!=0))
miss.values.per.class<- sapply(miss.per.class.usingRows,sum)
miss.df.per.class<- data.frame(class=names(miss.values.per.class),total=tot.values.per.class,missing=miss.values.per.class,missProp=miss.values.per.class/tot.values.per.class,stringsAsFactors=FALSE)
miss.df.per.class<- miss.df.per.class[order(miss.df.per.class$total),] miss.df.per.class$class<- factor(miss.df.per.class$class,levels=miss.df.per.class$class,ordered=FALSE)
melt.miss.df.class<- melt(miss.df.per.class[,c(1,2,3)],id.vars=1)
g3<- ggplot(data=melt.miss.df.class,aes(x=class,y=value)) + geom_bar(aes(fill=variable),position="identity") + labs(x="Classes",y="Values") + theme(axis.text.x=element_text(angle=45,hjust=1))


#now let's plot a heatmap for the proportion of missing values per classes against all variables

miss.per.class.usingCols<- sapply(df.per.class,function(x)apply(x,2,function(y)sum(is.na(y))))
miss.cols.per.class<- apply(miss.per.class.usingCols,2,function(y)sum(y!=0))
missing.melt<- melt(miss.per.class.usingCols)
values.per.class.usingCols<- sapply(df.per.class,function(x)apply(x,2,function(y)length(y)))
tot.values.melt<- melt(values.per.class.usingCols)
final.df<- data.frame(variables=missing.melt$Var1,classes=missing.melt$Var2,missing=missing.melt$value,total=tot.values.melt$value,missingProp=missing.melt$value/tot.values.melt$value)
g4<- ggplot(data=final.df,aes(x=classes,y=variables)) + geom_tile(aes(fill=missingProp)) + scale_fill_gradient(name="Proportion of missing values",low="white",high="red") + theme(axis.text.x=element_text(angle=45,hjust=1)) + labs(x="Classes",y="Variables")


From the above plot, I can see that there are certain variables that are missing for all 5 of the classes with missing data, but I am not wise as to whether these are MAR or MNAR.

I have used the following code to check this:

selected<- miss.per.class.usingCols[,c("2-4-d-injury","cyst-nematode","diaporthe-pod-&-stem-blight","herbicide-injury","phytophthora-rot")]
all.missing<- apply(selected,1,function(x)all(x!=0)==TRUE)
all.missing
#hail, sever, seed.tmt, germ, leaf.mild, lodging


The only step next that I could think of was to repeatedly generate mosaic plots for various combinations of variables to see if some relationship exists.

For example, I found that stem.cankers + canker.lesion, canker.lesion + ext.decay, and sever + seed.tmt seem to show some pattern in that all missing data for one variable occurs in cases of missing data for the other.

#inspect using mosaic plots for categorical data
d2<- data.frame(as.matrix(d1),stringsAsFactors=FALSE)
d2[is.na(d2)]<- "NAs"
d3<- data.frame(sapply(d2,as.factor))
mosaic(xtabs(~d3$stem.cankers + d3$canker.lesion))
mosaic(xtabs(~d3$sever + d3$seed.tmt))


But I don't find these plots useful in understanding whether my data is MAR or MNAR.

• Is there any better technique for my purpose that I am unaware of?
• Shall I keep working with mosaic plots to exhaustively check each pair of variables?
• How would you go about finding these relationships in data? Any pointers or advice would be appreciated.

I was not familiar with MAR (as opposed to MCAR i.e. Missing Completely At Random) and looked it up. Wikipedia's description is

Missing at random (MAR) occurs when the missingness is not random, but where missingness can be fully accounted for by variables where there is complete information.

The description continues, citing Roderick & Rubin (2002)

Since MAR is an assumption that is impossible to verify statistically, we must rely on its substantive reasonableness.

Little, Roderick J. A.; Rubin, Donald B. (2002), Statistical Analysis with Missing Data (2nd ed.), Wiley.

If this second claim is correct (and I have insufficient expertise to judge this), then trying to find a plot that will let you decide between MAR and MNAR is very unlikely to succeed.

For what it's worth, I would say that this plot of yours is clear evidence that the data is not MCAR: