R: Konstante in Linearkombination addieren, glht()

Question

Ich versuche also, eine Statafunktion zu replizieren, die ich in Principles of Econometrics sah, von Hill, Griffiths und Lim. Die Funktion, die ich replizieren möchte, sieht in Stata so aus; Diese

lincom _cons + b_1 * [arbitrary value] - c 

ist auf die Nullhypothese H0: B0 + B1 * X = C

Ich bin in der Lage, eine Hypothese ohne die Konstante zu testen, aber ich möchte die Konstante hinzufügen, wenn die lineare Kombination testen von Parametern. Ich habe die Paketdokumentation für glht() durchgesehen, aber es hatte nur ein Beispiel, wodurch sie die Konstante herausnahmen. Ich habe das Beispiel repliziert, wobei ich die Konstante beibehalten habe, aber ich bin unsicher, wie man eine lineare Kombination testen kann, wenn man eine Matrix K und eine Konstante hat. Als Referenz, hier ist ihr Beispiel;

### multiple linear model, swiss data 
lmod <- lm(Fertility ~ ., data = swiss) 

### test of H_0: all regression coefficients are zero 
### (ignore intercept) 

### define coefficients of linear function directly 
K <- diag(length(coef(lmod)))[-1,] 
rownames(K) <- names(coef(lmod))[-1] 
K 

### set up general linear hypothesis 
glht(lmod, linfct = K) 

Ich bin nicht wirklich gut in der Erstellung vorgetäuscht Datensätze, aber hier ist mein Versuch.

library(multcomp) 
test.data = data.frame(test.y = seq(200,20000,1000), 
        test.x = seq(10,1000,10)) 
test.data$test.y = sort(test.data$test.y + rnorm(100, mean = 10000, sd = 100)) - 
    rnorm(100, mean = 5733, sd = 77) 
test.lm = lm(test.y ~ test.x, data = test.data) 

# to view the name of the coefficients 
coef(test.lm) 

# this produces an error. How can I add this intercept? 
glht(test.lm, 
linfct = c("(Intercept) + test.x = 20")) 

Es scheint, dass es zwei Möglichkeiten gibt, um dies zu, die Dokumentation nach. Ich kann entweder die Funktion diag() verwenden, um eine Matrix zu erstellen, die ich dann im Argument linfct = verwenden kann, oder ich kann eine Zeichenfolge verwenden. Die Sache mit dieser Methode ist, dass ich nicht ganz weiß, wie man die Methode diag() benutzt, während man auch die Konstante (die rechte Seite der Gleichung) einbezieht; Im Fall der Zeichenkettenmethode bin ich nicht sicher, wie man den Achsenabschnitt hinzufügt.

Jede und alle Hilfe würde sehr geschätzt werden.

Hier sind die Daten, mit denen ich arbeite. Das war ursprünglich in einer .dta-Datei, also entschuldige ich mich für die schreckliche Formatierung. Laut dem oben erwähnten Buch ist dies die Datei food.dta.

structure(list(food_exp = structure(c(115.22, 135.98, 119.34, 
114.96, 187.05, 243.92, 267.43, 238.71, 295.94, 317.78, 216, 
240.35, 386.57, 261.53, 249.34, 309.87, 345.89, 165.54, 196.98, 
395.26, 406.34, 171.92, 303.23, 377.04, 194.35, 213.48, 293.87, 
259.61, 323.71, 275.02, 109.71, 359.19, 201.51, 460.36, 447.76, 
482.55, 438.29, 587.66, 257.95, 375.73), label = "household food expenditure per week", format.stata = "%10.0g"), 
income = structure(c(3.69, 4.39, 4.75, 6.03, 12.47, 12.98, 
14.2, 14.76, 15.32, 16.39, 17.35, 17.77, 17.93, 18.43, 18.55, 
18.8, 18.81, 19.04, 19.22, 19.93, 20.13, 20.33, 20.37, 20.43, 
21.45, 22.52, 22.55, 22.86, 24.2, 24.39, 24.42, 25.2, 25.5, 
26.61, 26.7, 27.14, 27.16, 28.62, 29.4, 33.4), label = "weekly household income", format.stata = "%10.0g")), .Names = c("food_exp","income"), class = c("tbl_df", "tbl", "data.frame"), row.names = c(NA, -40L)) 

Answer 1

Lassen Sie uns die Daten aus Ihrem Buch laden dann unsere Ergebnisse überprüfen, wie wir entlang gehen, um sicherzustellen, sind wir die gleiche Sache zu bekommen. Ich präsentiere Ihnen die Antwort teilweise, weil es mir geholfen hat, genau zu verstehen, was Sie wollten, und Sie teilweise von der Äquivalenz hier zu überzeugen.

Teil der Verwirrung für mich war in der Syntax Ihrer lincom Beispiel. Deine Syntax mag korrekt gewesen sein, ich habe keine Ahnung, aber basierend darauf, wie es aussah, dachte ich, du würdest etwas anderes machen, also war es hilfreich, auf dein Buch Bezug zu nehmen.

Lassen Sie uns zunächst die Daten laden und das lineare Modell ausführen, die auf Seite 115 ist:

install.packages("devtools") # if not already installed 
library(devtools) 
install_git("https://github.com/ccolonescu/PoEdata") 

library(PoEdata) # loads the package in memory 
library(multcomp) # for hypo testing 
data(food)   # loads the data set of interest 

# EDA 
summary(food) 

# Model 
mod <- lm(food_exp ~ income, data = food) 
summary(mod) # Note: same results as PoE 4th ed. Pg 115 (other than rounding) 

Call: 
lm(formula = food_exp ~ income, data = food) 

Residuals: 
    Min  1Q Median  3Q  Max 
-223.025 -50.816 -6.324 67.879 212.044 

Coefficients: 
      Estimate Std. Error t value Pr(>|t|)  
(Intercept) 83.416  43.410 1.922 0.0622 . 
income  10.210  2.093 4.877 1.95e-05 *** 
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 89.52 on 38 degrees of freedom 
Multiple R-squared: 0.385, Adjusted R-squared: 0.3688 
F-statistic: 23.79 on 1 and 38 DF, p-value: 1.946e-05 

So weit, so gut. Auf Pg. 115 der 4. Ausgabe zeigt das gleiche Regressionsmodell, abgesehen von einigen kleineren Rundungsdifferenzen.

Als nächstes wird das Buch eine Punktschätzung der wöchentlichen Lebensmittelausgaben berechnet, abhängig Haushaltseinkommen von 20 (die $ 2.000/Wo ist):

# Point estimate 
predict(mod, newdata = data.frame(income=20)) 

 1 
287.6089 

Wieder

, erhalten wir die genau dasselbe Ergebnis.Übrigens, Sie können auch die gleichen Ergebnisse in Stata in der schönen kostenlosen Probe des Buches Using Stata for Principles of Econometrics 4th ed. von Wiley sehen.

Endlich sind wir bereit für die Hypothesentests. Wie bereits erwähnt, möchte ich sicherstellen, dass ich genau reproduzieren kann, was Stata hatte. Sie haben Ihren Code freundlicherweise zur Verfügung gestellt, aber ich war etwas verwirrt über Ihre Syntax.

Zum Glück haben wir Glück gehabt. Während die Vorschau der 4.en Führung nur Stata Ausgabe 2 Kapitel durchläuft, waren die Economics and Business faculty eine Universität in Holland können Teile einer alten Ausgabe erhalten zur Verfügung gestellt DRM-frei, so können wir auf das verweisen:

siehe

und schließlich, dass wir es in R wie diese replizieren kann:

# Hyothesis Test 
summary(glht(mod, linfct = c("income = 15"))) 

Simultaneous Tests for General Linear Hypotheses 

Fit: lm(formula = food_exp ~ income, data = food) 

Linear Hypotheses: 
      Estimate Std. Error t value Pr(>|t|) 
income == 15 10.210  2.093 -2.288 0.0278 * 
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
(Adjusted p values reported -- single-step method) 

sie sich nicht von den verschie täuschen Mietausgabeformat. Die estimate zeigt Ihnen im R-Code ist einfach der Koeffizient ("b2") auf income aus dem Regressionsmodell. Sie ändert sich nicht in Abhängigkeit vom Hypothesentest, während sie in der Stata-Ausgabe "b2 - 15" (was in R mod$coefficients[2]-15 ist) tun.

Was sind die Werte t (t value) und p (Pr(>|t|))? Beachten Sie, dass die Teststatistiken von R denen von Stata entsprechen.

Es war ein weiteres Beispiel mit einem H0 Einkommen = 7,5 Mal sehen, dass der t-Wert 1,29 und p-Wert ist 0,203 sowohl R und Stata:

summary(glht(mod, linfct = c("income = 7.5"))) 

Simultaneous Tests for General Linear Hypotheses 

Fit: lm(formula = food_exp ~ income, data = food) 

Linear Hypotheses: 
       Estimate Std. Error t value Pr(>|t|) 
income == 7.5 10.210  2.093 1.294 0.203 
(Adjusted p values reported -- single-step method) 

Sie können die Konfidenzintervalle auch mit confint() abrufen.

Schließlich denke ich, Sie bei Abschnitt 3.6.4 (S. 117) Ihres Buches waren auf der Suche, wo ein Manager die Hypothese, dass, testen wollte income von 20 ($ 2000/wk) geben die food_exp ist> 250:

Wir können den T-Wert in R berechnen als:

t = sum((mod$coefficients[1] + 20*mod$coefficients[2]-250)/sqrt(vcov(mod)[1] + 20^2 * vcov(mod)[4] + 2 * 20 * vcov(mod)[2])) 
t 
[1] 2.652613 

Wo die Formel ist die gleiche wie auf den vorhergehenden Seiten 2 des Buchs.

Wir können auch diese in eine benutzerdefinierte Funktion machen (funktioniert für einfache lineare Regression, was bedeutet, 1 unabhängige Variable nur):

hypo_tester <- function(expenditure, income_per_week_hundreds, mod){ 
    t = sum((mod$coefficients[1] + 
      income_per_week_hundreds*mod$coefficients[2]-expenditure)/sqrt(vcov(mod)[1] + 
      income_per_week_hundreds^2 * vcov(mod)[4] + 2 * income_per_week_hundreds * vcov(mod)[2])) 
    return(t) 
} 

hypo_tester(250, 20, mod) 
[1] 2.652613 
hypo_tester(200, 20, mod) 
[1] 6.179193 
hypo_tester(300, 20, mod) 
[1] -0.8739668