I need some help with the statistical analysis of a study of a particular surgery to remove a particular cancer. I am using the statistical program R to conduct my analysis. My data are saved in the object study_data
.
Data
# Create reproducible example data
set.seed(50)
study_data <- data.frame(
Patient_ID = 1:500,
Institution = sample(c("New York","San Francisco","Houston","Chicago"),500,T),
Gender = sample(c("Male","Female"),500,T),
Race = sample(c("White","Black","Hispanic","Asian"),500,T),
Tumor_grade = sample(c("One","Two","Three","Four"),500,T),
Pathologic_stage = sample(c("P0","Pa","Pis","P1","P2a","P2b","P3a","P3b","P4a","P4b"),500,T),
Treatment_arm = sample(c("One","Two","Three","Four"),500,T),
Surgery_age = round(runif(500,20,100)),
Nodes_removed = round(runif(500,1,130)))
Here is what the data look like:
# Peak at the first six lines of the data
head(study_data)
Patient_ID Institution Gender Race Tumor_grade Pathologic_stage Treatment_arm Surgery_age Nodes_removed
1 1 Houston Male Hispanic One P2b Two 77 130
2 2 San Francisco Female Hispanic Three Pa Two 38 112
3 3 New York Female Black Four P0 Four 90 90
4 4 Chicago Male Hispanic Two Pis Four 46 4
5 5 Houston Female Black Four P2a Four 96 114
6 6 New York Male Black Three P3b Four 92 7
My interest
I am interested in learning more about what variables are associated with the number of lymph nodes removed during the surgery. My first thought was to simply stratify the data by a particular variable and then calculate the median number of nodes removed.
For example, to see if the institution at which the surgery was performed mattered, I could write:
cbind(do.call(rbind, by(study_data$Nodes_removed, study_data$Institution, summary)))
Min. 1st Qu. Median Mean 3rd Qu. Max.
Chicago 1 25.50 65.5 64.48 98.75 129
Houston 1 40.00 71.0 69.26 100.00 130
New York 4 36.00 67.0 67.96 100.00 129
San Francisco 3 36.75 61.0 65.76 99.00 127
This lets me compare the median nodes removed in each institutional city.
My question
I would like to fully examine the association between all of my variables and the outcome Nodes_removed
.
- Should I just do these simple summary statistics for all of my variables?
- Do I need to perform some sort of hypothesis test for all of the associations to say whether or not the summary statistics differ? For example, should I calculate a median and a confidence interval for each comparison?
- Or should I be using t-tests to compare one group to another?
- In the case of a multi-level variable, should I use ANOVA?
- Is there any role for linear regression analysis here?
- If I wanted to build a single model that includes every possible predictor variable, what method should I use?
For example, say that I am most interested in the association between the age at which the surgery was performed, Surgery_age
, and Nodes_removed
. However, I would like to adjust this association for potential confounders like gender, race, tumor grade, treatment arm, etc. What is the best way for me to do this?
Thanks for any advice you can give!
glm(..., family=poisson)
, b/c you have count data. I would start withsummary(glm(Nodes_removed~(.), data=study_data, family=poisson))
, but I notice the QQ plot doesn't look so great (if you plot the glm() object you get diagnostics). Also, there are a lot of terms in there, and I don't know what types of interactions etc. seem sensible. With data like this, I would do careful thinking about what makes sense, b/c blindly including a lot of terms can be dangerous. $\endgroup$