Um blog sobre nada

Um conjunto de inutilidades que podem vir a ser úteis

T – Test for regression Coefficients

Posted by Diego em Novembro 12, 2014


 

Getting the coefficients:

 

·         Prepare data:

library(UsingR);
data(diamond)

y <- diamond$price;
x <- diamond$carat;
n <- length(y)

 

·         Slope and intersect:

beta1 <- cor(y, x) * sd(y) / sd(x)

beta0 <- mean(y) – beta1 * mean(x)

·         Calculate the error:

e <- y – beta0 – beta1 * x #y – (beta1*x + beta0)

clip_image002[4]

·         Calculate Sigma

> sigma <- sqrt(sum(e^2) / (n-2))

> sigma

[1] 31.84052

·         Calculate the sum of the squared Xs:
> ssx <- sum((x – mean(x))^2)
>
ssx
[1] 0.1515667

·         Calculate the standard error for Beta0 and Beta1
> seBeta0 <- (1 / n + mean(x) ^ 2 / ssx) ^ .5 * sigma
>
seBeta1 <- sigma / sqrt(ssx)
>
seBeta0
[1] 17.31886
>
seBeta1
[1] 81.78588

·         T distribution for hypothesis testing:
> tBeta0 <- beta0 / seBeta0;
>
tBeta1 <- beta1 / seBeta1
>
tBeta0
[1] -14.99094
>
tBeta1
[1] 45.49715

·         P values (shouldn’t be bigger than 1):
> pBeta0 <- 2 * pt(abs(tBeta0), df = n – 2, lower.tail = FALSE)
>
pBeta1 <- 2 * pt(abs(tBeta1), df = n – 2, lower.tail = FALSE)
>
pBeta0
[1] 2.523271e-19
>
pBeta1
[1] 6.75126e-40

 

 

Similarly:

> fit<-lm(y~x);
> summary(fit)$coefficients

 

clip_image003[4]

 

Confidence Interval:

> sumCoef<-summary(fit)$coefficients
> sumCoef[2,1]+c(-1,1)*qt(.975,df=fit$df)*sumCoef[2,2]
[1] 3556.398 3885.651

 

With 95% confidence, we estimate that a 0.1 carat increase in diamond size results in a 355.6 to 388.6 increase in price in (Singapore) dollars.

In other words, the model predicts that a 0.16 carat diamond will cost 335.73, based on the confidence interval, we can assume with 95% confidence that a 0.26 (+0.1) diamond will cost between 691.33 (335.73 + 355.6) and 724.33 (335.73 + 388.6). This can be “verified” by looking at the predicted value for a 0.26 carat diamond, which is 707.83

 

Plotting the prediction Intervals:

Plot the points and the regression line:

> plot(x,y,frame=FALSE,xlab="Carat",ylab="Dollars",pch=21,col="black",bg="lightblue",cex=2)
> abline(fit,lwd=2)
 

 

Calculate X points and the predicted Y points:
 
> xVals<-seq(min(x),max(x),by=.01)
> yVals<-beta0+beta1*xVals

 

clip_image004[4]

 

Calculate the standard error:

> se1<-sigma*sqrt(1/n+(xVals-mean(x))^2/ssx)
> se2<-sigma*sqrt(1+1/n+(xVals-mean(x))^2/ssx

clip_image006[4]

 

> lines(xVals,yVals+2*se1)
> lines(xVals,yVals-2*se1)
> lines(xVals,yVals+2*se2)
> lines(xVals,yVals-2*se2)

 

 

clip_image008[4]

Deixe uma Resposta

Preencha os seus detalhes abaixo ou clique num ícone para iniciar sessão:

Logótipo da WordPress.com

Está a comentar usando a sua conta WordPress.com Terminar Sessão / Alterar )

Imagem do Twitter

Está a comentar usando a sua conta Twitter Terminar Sessão / Alterar )

Facebook photo

Está a comentar usando a sua conta Facebook Terminar Sessão / Alterar )

Google+ photo

Está a comentar usando a sua conta Google+ Terminar Sessão / Alterar )

Connecting to %s

 
%d bloggers like this: