**MeteorologistsBrier 2021-1021TN.dp
*Calculate Brier Score for Meterologist problem
*Tom Newman 10/21/21 

/* The first step is to get the data from the problem into a form Stata can understand.  You can enter the data in Excel and then import it, or enter it using Stata's data editor, or you can just the Stata dataset this do-file creates.  Since this is actually the hardest part, we might want to figure out how to enter the data yourself.  On the other hand usually your data will already be in a 1-record-per-observation form.*/
 
*If you import from Excel, you will have to change the path.

import excel "/Users/thomasnewman/Box/TEACH/CLINEPI/Clinical Epi 2021/Session 5 risk prediction/MeteorologistFrequencies.xls", sheet("Sheet1") firstrow clear

*I made the Excel table the same way as was illustrated when we were doing Kappa.  Here's what it looks like
bysort channel: list, clean
save meteorologists, replace

*Or you could skip the Excel importand just do this:
use meteorologists, clear

*We didn't do this next part when we did Kappa because Kappa accepts the [fweight=] option.  For Stata commands that don't, if you start with data in a table you first have to use the numbers in the COUNT field to EXPAND the dataset so that that each row represents one observation.  Usually these will be patients; in this example, they are days:

expand count /*This command duplicates each observation "count" times, to create a 1-record per day dataset for each channel*/
drop if count==0  /*This is needed because of how the expand command works*/
sort channel predict rain
*Look at the expanded dataset.  Note that now each observation represents a day for one of the two channels 
bysort channel: list channel predict rain

bysort channel: brier rain pred 

* Note: The "Sanders Resolution" is the same as the Refinement Loss.
* Note: "Reliability in the Small" is the same as the Calibration Error.

*Note that when I created my original spreadsheet, created the two channel predictions separately.  But if the data were only for 1 month, they are in LONG format, and it might be clearer to have them in WIDE format.  Following steps do this:

bysort channel:gen day=_n
list in 1/40
reshape wide prediction,i(day) j(channel)
list

*Now prediction2 and prediction3 are the predictions for those 2 channels:
brier rain prediction2
brier rain prediction3

*end