Kapitel 82 Stata in RMarkdown
Manchmal ist es hilfreich, Stata direkt in einem RMarkdown laufen zu lassen, zum Beispiel um Analysen zu vergleichen.
Stata muss natürlich auf dem Computer - mit einer gültigen Lizenz - installiert sein. Hier finden Sie eine Erklärung zur Installation.
82.0.1 Pakete für dieses Kapitel
Wir benötigen folgende Pakete:
library(Statamarkdown)
library(dplyr)82.1 Beispiel mit Graphik
- In RMarkdown fügen wir einen Stata code chunk ein (anstatt {r} steht da {stata})
Abbildung 82.1: Bild des Stata code chunks, mit dem wir die Graphik erstellen.
- Wir simulieren 500 Werte
- Stata zeichnet ein Histogramm
- Wir speichern dieses als .png
set obs 500
generate outcome = rnormal(100,12)
hist(outcome)
graph export "histogram_simulated_data.png", replace(bin=22, start=66.780296, width=3.0344561)
file histogram_simulated_data.png saved as PNG formatNun können wir dieses Histogramm in Rmarkdown anzeigen lassen:
knitr::include_graphics("histogram_simulated_data.png")
Abbildung 82.2: Graphik, die von Stata erstellt und gespeichert wurde, und danach von R als .png angezeigt wird.
82.2 Vergleich der Analysen R versus Stata
Wir laden in R Daten von einem Artikel herunten. Siehe hier für den Artikel:Guillaumier, A., Spratt, N. J., Pollack, M., Baker, A., Magin, P., Turner, A., … & Bonevski, B. (2022). Evaluation of an online intervention for improving stroke survivors’ health-related quality of life: A randomised controlled trial. PLoS Medicine, 19(4), e1003966..
Sie finden die benutzten Fragebogen hier (Baseline) und hier (Follow-up).
Wir wählen vier Variablen aus: ID, Raucherstatus bei Baseline, sowie eq5d Gesundheitsstatus bei Baseline und beim 6-Monats Follow-up.
Wir speichern die Daten im Stata .dta Format.
data<-rio::import("https://nova.newcastle.edu.au/vital/access/services/Download/uon:38922/ATTACHMENT01", format="xlsx")
data<-data %>%
select(ID, B_SMOKE,B_eq5d5lHealthState, eq5d5lHealthState )
rio::export(data, "data_Guillaumier_smoke_eq5d_health_status.dta")Nun können wir im RMarkdown die Daten in Stata laden und auswerten.
use "data_Guillaumier_smoke_eq5d_health_status.dta"
describe
save Guillaumier_from_Stata.dtaContains data from data_Guillaumier_smoke_eq5d_health_status.dta
Observations: 399
Variables: 4 24 Dec 2022 20:07
-------------------------------------------------------------------------------
Variable Storage Display Value
name type format label Variable label
-------------------------------------------------------------------------------
ID double %10.0g
B_SMOKE str36 %-9s
B_eq5d5lHealt~e str5 %-9s
eq5d5lHealthS~e str5 %-9s
-------------------------------------------------------------------------------
Sorted by:
file Guillaumier_from_Stata.dta already exists
r(602);
end of do-file
r(602);Das Problem ist, dass wir nicht einfach mehrere Code-Chunks machen können, da die Daten nicht übertragen werden. Wir müssen also die Daten in jedem Chunk neu laden. Zuerst müssen wir noch ein paar Variablen von String zu Numeric umwandeln.
tab B_SMOKE
destring B_eq5d5lHealthState eq5d5lHealthState , replace
encode B_SMOKE, gen(smoker)
tab B_SMOKE smokerr(111);
end of do-file
r(111);Also das gleiche noch einmal, diesmal mit laden der Daten:
use Guillaumier_from_Stata.dta
tab B_SMOKE
destring B_eq5d5lHealthState eq5d5lHealthState , replace
encode B_SMOKE, gen(smoker2)
tab B_SMOKE smoker2
save Guillaumier_from_Stata.dta, replace
describe B_SMOKE B_SMOKE | Freq. Percent Cum.
-------------------------------------+-----------------------------------
Not at all | 380 95.96 95.96
Yes, but less often than once a week | 1 0.25 96.21
Yes, daily | 15 3.79 100.00
-------------------------------------+-----------------------------------
Total | 396 100.00
B_eq5d5lHealthState already numeric; no replace
eq5d5lHealthState already numeric; no replace
variable smoker2 already defined
r(110);
end of do-file
r(110);Hier eine lineare Regression:
use Guillaumier_from_Stata.dta
regress eq5d5lHealthState i.smoker Source | SS df MS Number of obs = 355
-------------+---------------------------------- F(2, 352) = 0.47
Model | 52927147.5 2 26463573.7 Prob > F = 0.6255
Residual | 1.9825e+10 352 56322388 R-squared = 0.0027
-------------+---------------------------------- Adj R-squared = -0.0030
Total | 1.9878e+10 354 56153694.1 Root MSE = 7504.8
------------------------------------------------------------------------------
eq5d5lHeal~e | Coefficient Std. err. t P>|t| [95% conf. interval]
-------------+----------------------------------------------------------------
smoker |
Yes, but .. | -4893.849 7515.725 -0.65 0.515 -19675.22 9887.524
Yes, daily | 1714.651 2407.481 0.71 0.477 -3020.206 6449.508
|
_cons | 16004.85 404.633 39.55 0.000 15209.05 16800.65
------------------------------------------------------------------------------Jetzt können wir diese Analyse mit der Analyse in R vergleichen.
data<-rio::import("data_Guillaumier_smoke_eq5d_health_status.dta")
data<-data %>%
mutate(B_eq5d5lHealthState=as.numeric(B_eq5d5lHealthState),
eq5d5lHealthState =as.numeric(eq5d5lHealthState ),
smoker=as.factor(B_SMOKE))
table(data$smoker)
Not at all
3 380
Yes, but less often than once a week Yes, daily
1 15 options(width = 300)
lm.fit<-lm(eq5d5lHealthState~smoker, data=data)
summary(lm.fit)
Call:
lm(formula = eq5d5lHealthState ~ smoker, data = data)
Residuals:
Min 1Q Median 3Q Max
-6608 -4894 -4884 5206 27440
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 22211 7505 2.960 0.00329 **
smokerNot at all -6206 7516 -0.826 0.40950
smokerYes, but less often than once a week -11100 10613 -1.046 0.29635
smokerYes, daily -4492 7871 -0.571 0.56861
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 7505 on 352 degrees of freedom
(43 observations deleted due to missingness)
Multiple R-squared: 0.004565, Adjusted R-squared: -0.003919
F-statistic: 0.538 on 3 and 352 DF, p-value: 0.6565Wir beobachten hier etwas komisches. Wenn wir eine Faktorvariable haben mit drei Antwortsoptionen, dann zeigt R in der Regression nur 2 Variablen, die Referenz- oder Baselinevariable wird nicht gezeigt, genau gleich wie in Stata. Siehe das R Beispiel hier unten.
independent<-factor(sample(c("A", "B", "C"), size=500, rep=TRUE))
dependent<-rnorm(500, 50,12)
lm.fit<-lm(dependent~independent)
summary(lm.fit)
Call:
lm(formula = dependent ~ independent)
Residuals:
Min 1Q Median 3Q Max
-28.360 -8.636 -0.570 8.465 34.505
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 48.4726 0.9516 50.940 <2e-16 ***
independentB 0.1185 1.3414 0.088 0.930
independentC 0.1550 1.2853 0.121 0.904
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 11.85 on 497 degrees of freedom
Multiple R-squared: 3.098e-05, Adjusted R-squared: -0.003993
F-statistic: 0.0077 on 2 and 497 DF, p-value: 0.9923Im Beispiel zuvor, zeigte R jedoch alle drei Kategorien der Faktorvariable. Hmmm… komisch. Das Problem ist nun folgendes: Weil wir die Faktor-Variable etwas faul erstellt haben, ohne die Levels explizit anzugeben, hat R vier Kategorien erstellt. Das müssen wir nun ändern:
data<-rio::import("data_Guillaumier_smoke_eq5d_health_status.dta")
data<-data %>%
mutate(B_eq5d5lHealthState=as.numeric(B_eq5d5lHealthState),
eq5d5lHealthState =as.numeric(eq5d5lHealthState ),
smoker=factor(B_SMOKE, levels=c("Not at all", "Yes, but less often than once a week", "Yes, daily")))
table(data$smoker)
Not at all Yes, but less often than once a week Yes, daily
380 1 15 So, jetzt sieht es besser aus. Schauen wir uns also noch einmal die lineare Regresson an:
options(width = 300)
lm.fit<-lm(eq5d5lHealthState~smoker, data=data)
summary(lm.fit)
Call:
lm(formula = eq5d5lHealthState ~ smoker, data = data)
Residuals:
Min 1Q Median 3Q Max
-6608 -4894 -4884 5206 27440
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 16004.8 404.6 39.554 <2e-16 ***
smokerYes, but less often than once a week -4893.8 7515.7 -0.651 0.515
smokerYes, daily 1714.7 2407.5 0.712 0.477
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 7505 on 352 degrees of freedom
(44 observations deleted due to missingness)
Multiple R-squared: 0.002663, Adjusted R-squared: -0.003004
F-statistic: 0.4699 on 2 and 352 DF, p-value: 0.6255Wir sehen hier, dass die Resultate nun vergleichbar mit Stata sind.