What tool should I use to find improper relationship between doctor and patient? I have a data table in R where I summarized the dollar value of prescribed medicines for each doctor-patient combination. In all, I have about half a million observations, with about 3000 doctors and 250 thousand patients (3 column data table: DoctorID, PatientID, totalValue).
I would like to find, based on data, if there is any indication of any of the doctors prescribing unusually (from the perspective of their own prescription habits) high value formulas to any specific patient (most likely any given doctor would have prescribed drugs to more than one patient. Similarily, any given patient would have been prescribed by more than one doctor), and thus find "improper" associations.
I thought ANOVA could be a tool to perform this analysis, but I don't know if this assumption is right, and if it's so, I don't know how to interpret the results (I ran aov(totalValue ~ DoctorID * PatientID,DataTable) and got Pr>F smaller than 1E-17 for both variables, and about 0.01 for the interaction, yet I don't have a clue as to how I can use those results to find Doctor-Patient pairs who could be misbehaving).
I appreciate if you can guide me as to what tool to use, and then how to interpret the produced results.
I'll also appreciate if someone more experienced could help me adding appropriate tags to this question. I'm so lost I can't even think of adequate tags.
============== EDIT TO ADD ===============
Please find a scaled-down, modified version of my data here. Each observation is one medicine from a formula.
I can obtain the sum of Value grouped by FormulaNumber for each formula, along with the prescribing Doctor and the Patient:
ex_doc_pat<-exDisp[,.(DoctorID=unique(DoctorID),PatientID=unique(PatientID),Disease=unique(Disease),totalValue=sum(Value)),by=FormulaNumber]

I can produce a 2way table from there (which, as @kjetil metioned, is fairly sparse --actual data is less sparse than this scaled-down version--):
2wayTab<-reshape(ex_doc_pat[,.(DoctorID,PatientID,totalValue)],idvar = "DoctorID",timevar = "PatientID",direction="wide")

Responding to @Scortchi, I dropped the rest of data, as it wasn't particularily interesting (place of the farmacy, code of the farmacy, etc); and yes, it could be that a high-tag prescription is well deserved by patient's condition. I just need a starting place to begin looking for unusually high priced doctor-patient relationships.
 A: Interesting question!
As a commenter pointed out, you will not be able to prove anything using these data as there may be other explanations, but I think you can use the data to identify a group of suspect doctors.
I assume that this is about prescription sedatives/analgesics or something similar so that you want to identify doctors that prescribe high doses to addicts.
I would start by defining the "improper relationship", whatever that means? I would mostly be interested in trying to identify those doctors who consistently prescribes high doses of these drugs. A high dose might be defined as being more than the 95:th percentile of what patients are prescribed, so for each doctor-patient-drug combination (should be only one of each, right?) we could create a new binary variable that is 0 if the dose (Value) is lower than the 95th percentile for that drug and 1 if higher. Now we could count the frequency of 1:s on this variable for each doctor, and we should find that some doctors will have very few, and some will have many. We could call this percentage variable "ProportionHigh" to reflect that it is the proportion of "too high" doses prescribed, and we give this proportion to each doctor
Now we would like to identify the doctors that have a high "ProportionHigh" value, but the proportion might be influenced by Speciality and Disease, both of which could be included as random effects in the model.
So we might get a mixed model with something like: ProportionHigh ~ 1 + (1|Speciality) + (1|Disease). I don't think it makes sense to add Drug as a random effect, as per definition, all drugs will have a 5% probability of having a 1 and 95% probability of having a 0. Another problem is what distribution function to use. Most doctors will have a very low ProportionHigh value, so a linear model is clearly not an option. Perhaps a proportional binomial model can be fitted, but I fear this might violate the assumption of independence between the observed values for each doctor. This would assume that if a doctor sees 30 patients, the prescription to one patient does not affect any prescription to other patients - so this assumes that each doctor only prescribes one drug to each patient. Perhaps an idea is to run separate analyses for separate drugs then?
EDIT: I just saw that there are 245 drugs so separate analyses per drug is out of the question. But perhaps this could be dealt with by including patient as a random effect?
Anyway, when we have our regression model, we can calculate the residuals for each row in the data. A high residual should mean that the row (doctor/patient/drug combination) has a higher score than predicted from the model. You could then calculate the proportion of high residuals for each doctor, say those over the 95th percentile again, and you will then have a list of the 5% of doctors who prescribe the highest doses when taking the speciality and disease into account.
All of this is really just me shooting from the hip, and I'm sure my method lacks statistical rigor. The 95th percentiles are of course arbitrary, and a high dose might be defined otherwise (there are usually maximum recommended doses for each drug, so a higher dose than this might be another way to define "too high").
A: @PavoDive, you have a pretty interesting data. Usually I like very much visualize the data and after studding correlations. So, maybe is not worth doing but, I would create 2 new sets:
- patientId; number of doctors he/she visit; total of amount associated with the patient
- doctorId; number of patients he/she prescribed; total amount of doctors prescriptions
After I would do a scatter plot and check for anything "different". You can also start for correlating data.
After that I would start "digging": behavior of patients with more than one doctor, you can calculate the mean value for both total doctors/patients prescriptions, check who is above the mean, etc... 
Just start by knowing your data before you get into more complex analysis. 
A: ##### NOTE TO READERS #####
This is what I've done this far, so I'm posting it as an answer hoping that some "illuminated" will come and drop comments that allow me to learn more about the subject. That said, please comment and signal wherever you think there's a blatant error, whether conceptual or procedimental.



*

*From the data I built a data frame with the following fields: diagnostic (disease), doctor, patient and value for all doctor-patient pairs with existing values.

*For each Diagnostic I assumed the values of formulas (from any doctor) follow a normal distribution, so I can assign each case a probability depending on where it sits in comparison with the other values for the same disease.

*I estimated the geometric mean of those case-probabilities for each doctor-patient pair (say, for example: doctor A prescribed patient B $20 to treat diesase C, and given that it's an outlier for that specific disease, the probability is 0.01. However, A prescribed B a formula with probability 0.9 for disease D. The geomtric mean would be (0.01 x 0.9)^(1/2) = 0.095

*Low geometric means (specially those for large number of formulas) would suggest that that specific doctor usually prescribe drugs to that specific patient with values that are outliers to the distribution of values for the same diagnostic, thus rising a warning.


In code (with the data I provided in the question):
library(data.table)    

newDisp<-unique(exDisp[,
  .(concat=paste(DoctorID,PatientID,sep=""),
    Sumvalue=sum(Value),Disease=Disease),
  by=.(FormulaNumber)])

newDisp[,
   prob:=pnorm(Sumvalue,mean(Sumvalue),sd(Sumvalue),lower.tail = FALSE),
   by=Disease]

geomean<-function(x){prod(x)^(1/length(x))}

newDisp[,
  .(cumprob=geomean(prob),numForm=length(prob)),
  by=concat]
  [numForm>=3]
  [order(cumprob)]

which, for the toy data I provided doesn't show any interesting results (actual data has more than 2 million formulas, so it's likely to have significantly higher numForm values), but it does for the real data (very low cumprob for relatively high numForm values)
Finding these small cumprobs doesn't automatically mean it's fishy, but certainly that it smells so, so digging further is a good idea.
Of course I see there are a lot of problems with this approach. First, the normality assumption won't hold fast for many of the diseases, specially for those rare, "orphan" diseases where there could be too few points to even think of a distribution. Then, the geometric mean could be strongly affected by one value (which may not even be "improper"), thus signaling some false-positives. Finally, the data could be fooled by the system, say if a doctor is one of a few in her specialty and does wrong (prescribes high) with all his patients, thus creating a higher-than-real distribution that, nonetheless, would still be normal (this would be rather unlikely, as there would be too many people to involve in a criminal act). I appreciate if you can point more weaknesses of this approach.
