Styjun

In the movie Harry Potter and the Order or the Phoenix, why didn't Mr. Filch succeed to open the Room of Requirement if it's what he needed?

Why couldn't soldiers sight their own weapons without officers' orders?

Decode a variable-length quantity

How to help new students accept function notation

Could one become a successful researcher by writing some really good papers while being outside academia?

How to avoid ci-driven development..?

How can I tell if a flight itinerary is fake

How do I get the =LEFT function in excel, to also take the number zero as the first number?

If there were no space agencies, could a person go to space?

Why should public servants be apolitical?

4-dimensional Knight's Tour

Does this smartphone photo show Mars just below the Sun?

How to realistically deal with a shield user?

How does The Fools Guild make its money?

How would a family travel from Indiana to Texas in 1911?

What is the German idiom or expression for when someone is being hypocritical against their own teachings?

How can glass marbles naturally occur in a desert?

Colleagues speaking another language and it impacts work

one fifteen euros=115 euros: Is there a german equivalent?

Casting Goblin Matron with Plague Engineer on the battlefield

What was the first multiprocessor x86 motherboard?

Double blind peer review when paper cites author's GitHub repo for code

If a pawn is promoted to queen beside the king is the king instantly in check?

Do the Kraus operators of a CPTP channel need to be orthogonal?

Calculating Confidence intervals in R

adjusted bootstrap confidence intervals (BCa) with parametric bootstrap in boot packageConfidence intervals from lme's “intervals” on log plot (R)Stack sets of three columns with means and confidence intervals for multiple samplesReplicate least squares regression to check consistency of estimation and prediction with truthConfidence interval of polynomial regressionStargazer Confidence Interval Incorrect?Manually calculating the confidence interval of a multiple linear regression(OLS)R: Confidence intervals on non-linear fit with a non-analytic modelCalculating confidence Interval around a Linear Regression LineCalculate 95% confidence interval on the mean

.everyoneloves__top-leaderboard:empty,.everyoneloves__mid-leaderboard:empty,.everyoneloves__bot-mid-leaderboard:empty margin-bottom:0;

2

I'm trying to calculate confidence intervals from a t test in R manually and I suspect the way i calculate them are off.

Here's how I calculate confidence intervals manually right now

library(broom)
data("mtcars")
a1=tidy(t.test(mpg ~ am, mtcars))
mean_diff<-a1$estimate 
tvalue <-a1$statistic

#standard error 
sd1=sd(mtcars$mpg[mtcars$am==0])
sd2=sd(mtcars$mpg[mtcars$am==1])
n1=length(mtcars$mpg[mtcars$am==0])
n2=length(mtcars$mpg[mtcars$am==1])
#formula for standard error
stan_error=sqrt((sd1/n1)+(sd2/n2))

then I take the formula from this page about calculating the confidence intervals http://onlinestatbook.com/2/estimation/difference_means.html

The lower confidence interval I calculate like this

lower=mean_diff - (tvalue * stan_error)'

and the result comes out to be -4.147333

But the confidence intervals of

t.test(mpg ~ am, mtcars)

are

95 percent confidence interval:
 -11.280194 -3.209684

Any ideas?

r

share|improve this question

asked Mar 27 at 5:55

JLRJLR

4044 bronze badges

add a comment | 

2

I'm trying to calculate confidence intervals from a t test in R manually and I suspect the way i calculate them are off.

Here's how I calculate confidence intervals manually right now

library(broom)
data("mtcars")
a1=tidy(t.test(mpg ~ am, mtcars))
mean_diff<-a1$estimate 
tvalue <-a1$statistic

#standard error 
sd1=sd(mtcars$mpg[mtcars$am==0])
sd2=sd(mtcars$mpg[mtcars$am==1])
n1=length(mtcars$mpg[mtcars$am==0])
n2=length(mtcars$mpg[mtcars$am==1])
#formula for standard error
stan_error=sqrt((sd1/n1)+(sd2/n2))

then I take the formula from this page about calculating the confidence intervals http://onlinestatbook.com/2/estimation/difference_means.html

The lower confidence interval I calculate like this

lower=mean_diff - (tvalue * stan_error)'

and the result comes out to be -4.147333

But the confidence intervals of

t.test(mpg ~ am, mtcars)

are

95 percent confidence interval:
 -11.280194 -3.209684

Any ideas?

r

share|improve this question

asked Mar 27 at 5:55

JLRJLR

4044 bronze badges

add a comment | 

I'm trying to calculate confidence intervals from a t test in R manually and I suspect the way i calculate them are off.

Here's how I calculate confidence intervals manually right now

library(broom)
data("mtcars")
a1=tidy(t.test(mpg ~ am, mtcars))
mean_diff<-a1$estimate 
tvalue <-a1$statistic

#standard error 
sd1=sd(mtcars$mpg[mtcars$am==0])
sd2=sd(mtcars$mpg[mtcars$am==1])
n1=length(mtcars$mpg[mtcars$am==0])
n2=length(mtcars$mpg[mtcars$am==1])
#formula for standard error
stan_error=sqrt((sd1/n1)+(sd2/n2))

then I take the formula from this page about calculating the confidence intervals http://onlinestatbook.com/2/estimation/difference_means.html

The lower confidence interval I calculate like this

lower=mean_diff - (tvalue * stan_error)'

and the result comes out to be -4.147333

But the confidence intervals of

t.test(mpg ~ am, mtcars)

are

95 percent confidence interval:
 -11.280194 -3.209684

Any ideas?

r

share|improve this question

asked Mar 27 at 5:55

JLRJLR

4044 bronze badges

library(broom)
data("mtcars")
a1=tidy(t.test(mpg ~ am, mtcars))
mean_diff<-a1$estimate 
tvalue <-a1$statistic

#standard error 
sd1=sd(mtcars$mpg[mtcars$am==0])
sd2=sd(mtcars$mpg[mtcars$am==1])
n1=length(mtcars$mpg[mtcars$am==0])
n2=length(mtcars$mpg[mtcars$am==1])
#formula for standard error
stan_error=sqrt((sd1/n1)+(sd2/n2))

lower=mean_diff - (tvalue * stan_error)'

t.test(mpg ~ am, mtcars)

95 percent confidence interval:
 -11.280194 -3.209684

share|improve this question

asked Mar 27 at 5:55

JLRJLR

4044 bronze badges

share|improve this question

asked Mar 27 at 5:55

JLRJLR

4044 bronze badges

asked Mar 27 at 5:55

JLRJLR

4044 bronze badges

asked Mar 27 at 5:55

JLRJLR

4044 bronze badges

1 Answer
1

active

oldest

votes

1

Firstly, critical value for t is not right.

tvalue <- a1$statistic needs to be replaced with tvalue <- abs(qt(0.05/2, 30)).

Note its not 32 because we lose 2 degrees of freedom.

And you are missing ^2 (i.e. to the power of two) in formula for standard error.
What you have in sd1 and sd2 are standard errors so you need to convert this into variances.
So correct formula is:

stan_error = sqrt((sd1^2 / n1) + (sd2^2 / n2))

So the new code becomes:

library(broom)
data("mtcars")
a1=tidy(t.test(mpg ~ am, mtcars))
mean_diff<-a1$estimate 

t_cv<- abs(qt(0.05/2, 30))

#standard error 
sd1=sd(mtcars$mpg[mtcars$am==0])
sd2=sd(mtcars$mpg[mtcars$am==1])
n1=length(mtcars$mpg[mtcars$am==0])
n2=length(mtcars$mpg[mtcars$am==1])

#formula for standard error
stan_error = sqrt((sd1^2 / n1) + (sd2^2 / n2))

lower=mean_diff - (t_cv* stan_error)

lower
[1] -11.17264

But this still doesn't match the confidence interval using t.test function because t.test uses Welch's t-test (https://en.wikipedia.org/wiki/Welch%27s_t-test)
So your t-critical value under Welch's t-test should be

# Welch's t test degrees of freedom
welch_df <- (sd1^2/n1 + sd2^2/n2)^2 / (sd1^4/(n1^2*(n1-1)) + sd2^4/(n2^2*(n2-1)))
t_cv <- abs(qt(0.05/2, welch_df)) 

# Recalculate lower confidence interval
lower= mean_diff - (t_cv* stan_error)
lower
[1] -11.28019 # this matches confidence interval in t.test

share|improve this answer

edited Mar 27 at 7:12

answered Mar 27 at 6:55

MKaMKa

84677 silver badges1616 bronze badges

• 
  
  
  

  
  1
  

  
  
  
  
  
  

  
  
  So what even is this t value or "t" calculated by t.test (I got -3.7671) if it's different than the wtvalue (I got 18.33225 )and tvalue (I got 2.098196) which are actually using the real welch t test formulas to calculate the t value?
  
  – JLR
  Mar 27 at 7:03
  
  

  
  
  
  


• 
  
  
  

  
  1
  

  
  
  
  
  
  

  
  
  t statistic (or t value) is different to t critical value. t critical value is for significance testing, i.e. benchmark. t value is your calculated test statistic to compare against the critical value. You can get t value by mean_diff/stan_error = 3.7671
  
  – MKa
  Mar 27 at 7:09
  
  

  
  
  
  


• 
  
  
  

  
  1
  

  
  
  
  
  
  

  
  
  I think it is the object name that is confusing, just renamed tvalue to t_cv (t critical value)
  
  – MKa
  Mar 27 at 7:13
  
  

  
  
  
  

add a comment | 

Your Answer

StackExchange.ifUsing("editor", function ()
StackExchange.using("externalEditor", function ()
StackExchange.using("snippets", function ()
StackExchange.snippets.init();
);
);
, "code-snippets");

StackExchange.ready(function()
var channelOptions =
tags: "".split(" "),
id: "1"
;
initTagRenderer("".split(" "), "".split(" "), channelOptions);

StackExchange.using("externalEditor", function()
// Have to fire editor after snippets, if snippets enabled
if (StackExchange.settings.snippets.snippetsEnabled)
StackExchange.using("snippets", function()
createEditor();
);

else
createEditor();

);

function createEditor()
StackExchange.prepareEditor(
heartbeatType: 'answer',
autoActivateHeartbeat: false,
convertImagesToLinks: true,
noModals: true,
showLowRepImageUploadWarning: true,
reputationToPostImages: 10,
bindNavPrevention: true,
postfix: "",
imageUploader:
brandingHtml: "Powered by u003ca class="icon-imgur-white" href="https://imgur.com/"u003eu003c/au003e",
contentPolicyHtml: "User contributions licensed under u003ca href="https://creativecommons.org/licenses/by-sa/3.0/"u003ecc by-sa 3.0 with attribution requiredu003c/au003e u003ca href="https://stackoverflow.com/legal/content-policy"u003e(content policy)u003c/au003e",
allowUrls: true
,
onDemand: true,
discardSelector: ".discard-answer"
,immediatelyShowMarkdownHelp:true
);



);

draft saved

draft discarded

Sign up or log in

StackExchange.ready(function ()
StackExchange.helpers.onClickDraftSave('#login-link');
);

Sign up using Google

Sign up using Facebook

Sign up using Email and Password

Post as a guest

Name

Email

Required, but never shown

StackExchange.ready(
function ()
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fstackoverflow.com%2fquestions%2f55370644%2fcalculating-confidence-intervals-in-r%23new-answer', 'question_page');

);

Post as a guest

Name

Email

Required, but never shown

1

Firstly, critical value for t is not right.

tvalue <- a1$statistic needs to be replaced with tvalue <- abs(qt(0.05/2, 30)).

Note its not 32 because we lose 2 degrees of freedom.

And you are missing ^2 (i.e. to the power of two) in formula for standard error.
What you have in sd1 and sd2 are standard errors so you need to convert this into variances.
So correct formula is:

stan_error = sqrt((sd1^2 / n1) + (sd2^2 / n2))

So the new code becomes:

library(broom)
data("mtcars")
a1=tidy(t.test(mpg ~ am, mtcars))
mean_diff<-a1$estimate 

t_cv<- abs(qt(0.05/2, 30))

#standard error 
sd1=sd(mtcars$mpg[mtcars$am==0])
sd2=sd(mtcars$mpg[mtcars$am==1])
n1=length(mtcars$mpg[mtcars$am==0])
n2=length(mtcars$mpg[mtcars$am==1])

#formula for standard error
stan_error = sqrt((sd1^2 / n1) + (sd2^2 / n2))

lower=mean_diff - (t_cv* stan_error)

lower
[1] -11.17264

But this still doesn't match the confidence interval using t.test function because t.test uses Welch's t-test (https://en.wikipedia.org/wiki/Welch%27s_t-test)
So your t-critical value under Welch's t-test should be

# Welch's t test degrees of freedom
welch_df <- (sd1^2/n1 + sd2^2/n2)^2 / (sd1^4/(n1^2*(n1-1)) + sd2^4/(n2^2*(n2-1)))
t_cv <- abs(qt(0.05/2, welch_df)) 

# Recalculate lower confidence interval
lower= mean_diff - (t_cv* stan_error)
lower
[1] -11.28019 # this matches confidence interval in t.test

share|improve this answer

edited Mar 27 at 7:12

answered Mar 27 at 6:55

MKaMKa

84677 silver badges1616 bronze badges

• 
  
  
  

  
  1
  

  
  
  
  
  
  

  
  
  So what even is this t value or "t" calculated by t.test (I got -3.7671) if it's different than the wtvalue (I got 18.33225 )and tvalue (I got 2.098196) which are actually using the real welch t test formulas to calculate the t value?
  
  – JLR
  Mar 27 at 7:03
  
  

  
  
  
  


• 
  
  
  

  
  1
  

  
  
  
  
  
  

  
  
  t statistic (or t value) is different to t critical value. t critical value is for significance testing, i.e. benchmark. t value is your calculated test statistic to compare against the critical value. You can get t value by mean_diff/stan_error = 3.7671
  
  – MKa
  Mar 27 at 7:09
  
  

  
  
  
  


• 
  
  
  

  
  1
  

  
  
  
  
  
  

  
  
  I think it is the object name that is confusing, just renamed tvalue to t_cv (t critical value)
  
  – MKa
  Mar 27 at 7:13
  
  

  
  
  
  

add a comment | 

1

Firstly, critical value for t is not right.

tvalue <- a1$statistic needs to be replaced with tvalue <- abs(qt(0.05/2, 30)).

Note its not 32 because we lose 2 degrees of freedom.

And you are missing ^2 (i.e. to the power of two) in formula for standard error.
What you have in sd1 and sd2 are standard errors so you need to convert this into variances.
So correct formula is:

stan_error = sqrt((sd1^2 / n1) + (sd2^2 / n2))

So the new code becomes:

library(broom)
data("mtcars")
a1=tidy(t.test(mpg ~ am, mtcars))
mean_diff<-a1$estimate 

t_cv<- abs(qt(0.05/2, 30))

#standard error 
sd1=sd(mtcars$mpg[mtcars$am==0])
sd2=sd(mtcars$mpg[mtcars$am==1])
n1=length(mtcars$mpg[mtcars$am==0])
n2=length(mtcars$mpg[mtcars$am==1])

#formula for standard error
stan_error = sqrt((sd1^2 / n1) + (sd2^2 / n2))

lower=mean_diff - (t_cv* stan_error)

lower
[1] -11.17264

But this still doesn't match the confidence interval using t.test function because t.test uses Welch's t-test (https://en.wikipedia.org/wiki/Welch%27s_t-test)
So your t-critical value under Welch's t-test should be

# Welch's t test degrees of freedom
welch_df <- (sd1^2/n1 + sd2^2/n2)^2 / (sd1^4/(n1^2*(n1-1)) + sd2^4/(n2^2*(n2-1)))
t_cv <- abs(qt(0.05/2, welch_df)) 

# Recalculate lower confidence interval
lower= mean_diff - (t_cv* stan_error)
lower
[1] -11.28019 # this matches confidence interval in t.test

share|improve this answer

edited Mar 27 at 7:12

answered Mar 27 at 6:55

MKaMKa

84677 silver badges1616 bronze badges

• 
  
  
  

  
  1
  

  
  
  
  
  
  

  
  
  So what even is this t value or "t" calculated by t.test (I got -3.7671) if it's different than the wtvalue (I got 18.33225 )and tvalue (I got 2.098196) which are actually using the real welch t test formulas to calculate the t value?
  
  – JLR
  Mar 27 at 7:03
  
  

  
  
  
  


• 
  
  
  

  
  1
  

  
  
  
  
  
  

  
  
  t statistic (or t value) is different to t critical value. t critical value is for significance testing, i.e. benchmark. t value is your calculated test statistic to compare against the critical value. You can get t value by mean_diff/stan_error = 3.7671
  
  – MKa
  Mar 27 at 7:09
  
  

  
  
  
  


• 
  
  
  

  
  1
  

  
  
  
  
  
  

  
  
  I think it is the object name that is confusing, just renamed tvalue to t_cv (t critical value)
  
  – MKa
  Mar 27 at 7:13
  
  

  
  
  
  

add a comment | 

Firstly, critical value for t is not right.

tvalue <- a1$statistic needs to be replaced with tvalue <- abs(qt(0.05/2, 30)).

Note its not 32 because we lose 2 degrees of freedom.

And you are missing ^2 (i.e. to the power of two) in formula for standard error.
What you have in sd1 and sd2 are standard errors so you need to convert this into variances.
So correct formula is:

stan_error = sqrt((sd1^2 / n1) + (sd2^2 / n2))

So the new code becomes:

library(broom)
data("mtcars")
a1=tidy(t.test(mpg ~ am, mtcars))
mean_diff<-a1$estimate 

t_cv<- abs(qt(0.05/2, 30))

#standard error 
sd1=sd(mtcars$mpg[mtcars$am==0])
sd2=sd(mtcars$mpg[mtcars$am==1])
n1=length(mtcars$mpg[mtcars$am==0])
n2=length(mtcars$mpg[mtcars$am==1])

#formula for standard error
stan_error = sqrt((sd1^2 / n1) + (sd2^2 / n2))

lower=mean_diff - (t_cv* stan_error)

lower
[1] -11.17264

But this still doesn't match the confidence interval using t.test function because t.test uses Welch's t-test (https://en.wikipedia.org/wiki/Welch%27s_t-test)
So your t-critical value under Welch's t-test should be

# Welch's t test degrees of freedom
welch_df <- (sd1^2/n1 + sd2^2/n2)^2 / (sd1^4/(n1^2*(n1-1)) + sd2^4/(n2^2*(n2-1)))
t_cv <- abs(qt(0.05/2, welch_df)) 

# Recalculate lower confidence interval
lower= mean_diff - (t_cv* stan_error)
lower
[1] -11.28019 # this matches confidence interval in t.test

share|improve this answer

edited Mar 27 at 7:12

answered Mar 27 at 6:55

MKaMKa

84677 silver badges1616 bronze badges

stan_error = sqrt((sd1^2 / n1) + (sd2^2 / n2))

library(broom)
data("mtcars")
a1=tidy(t.test(mpg ~ am, mtcars))
mean_diff<-a1$estimate 

t_cv<- abs(qt(0.05/2, 30))

#standard error 
sd1=sd(mtcars$mpg[mtcars$am==0])
sd2=sd(mtcars$mpg[mtcars$am==1])
n1=length(mtcars$mpg[mtcars$am==0])
n2=length(mtcars$mpg[mtcars$am==1])

#formula for standard error
stan_error = sqrt((sd1^2 / n1) + (sd2^2 / n2))

lower=mean_diff - (t_cv* stan_error)

lower
[1] -11.17264

# Welch's t test degrees of freedom
welch_df <- (sd1^2/n1 + sd2^2/n2)^2 / (sd1^4/(n1^2*(n1-1)) + sd2^4/(n2^2*(n2-1)))
t_cv <- abs(qt(0.05/2, welch_df)) 

# Recalculate lower confidence interval
lower= mean_diff - (t_cv* stan_error)
lower
[1] -11.28019 # this matches confidence interval in t.test

share|improve this answer

edited Mar 27 at 7:12

answered Mar 27 at 6:55

MKaMKa

84677 silver badges1616 bronze badges

answered Mar 27 at 6:55

MKaMKa

84677 silver badges1616 bronze badges

answered Mar 27 at 6:55

MKaMKa

84677 silver badges1616 bronze badges