options nodate nocenter ps=1000; * ----------------------------------------------; * Exemplo 8.6.2 da 1a edição do livro do Rencher; * ----------------------------------------------; data ChemReaction; * x1=temperature x2=concentration x3=time y1=unchanged y2=converted; input x1 x2 x3 y1 y2; cards; 162 23.0 3.0 41.5 45.9 162 23.0 8.0 33.8 53.3 162 30.0 5.0 27.7 57.5 162 30.0 8.0 21.7 58.8 172 25.0 5.0 19.9 60.6 172 25.0 8.0 15.0 58.0 172 30.0 5.0 12.2 58.6 172 30.0 8.0 4.3 52.4 167 27.5 6.5 19.3 56.9 177 27.5 6.5 6.4 55.4 157 27.5 6.5 37.6 46.9 167 32.5 6.5 18.0 57.3 167 22.5 6.5 26.3 55.0 167 27.5 9.5 9.9 58.9 167 27.5 3.5 25.0 50.3 177 20.0 6.5 14.1 61.1 177 20.0 6.5 15.2 62.9 160 34.0 7.5 15.9 60.0 160 34.0 7.5 19.6 60.6 ; proc iml; *pág.292; use ChemReaction; read all var{x1} into x1; read all var{x2} into x2; read all var{x3} into x3; read all var{y1} into y1; read all var{y2} into y2; y = y2; * Exemplo 8.6.2; n = nrow(y); jn = j(n,1,1); In = I(n); X = jn||x1||x2||x3; k = 3; * número de variáveis regressoras; Beta = inv(t(X)*X)*t(X)*y; /* C = {0 1 -1 0, 0 0 2 -1}; gl_H0 = nrow(C); CBeta = C*Beta; SQH0 = t(CBeta)*inv(C*(inv(t(X)*X))*t(C))*CBeta; QMH0 = SQH0/gl_H0; SQRes = t(y)*(In-X*inv(t(X)*X)*t(X))*y; gl_res = n-k-1; QMRes = SQRes/gl_res; */ y = y1; * <<<<<< Seleciona a variável a ser analisada; Beta = inv(t(X)*X)*t(X)*y; G = inv(t(X)*X); a0={1,0,0,0}; Beta0 = t(a0)*Beta; G00 = t(a0)*G*a0; a1={0,1,0,0}; Beta1 = t(a1)*Beta; G11 = t(a1)*G*a1; a2={0,0,1,0}; Beta2 = t(a2)*Beta; G22 = t(a2)*G*a2; a3={0,0,0,1}; Beta3 = t(a3)*Beta; G33 = t(a3)*G*a3; SQRes = t(y)*(In-X*inv(t(X)*X)*t(X))*y; gl_res = n-k-1; QMRes = SQRes/gl_res; s = sqrt(QMRes); t_tab = tinv(0.975,n-k-1); * Calcula t-tabelado para alfa = 0,05/2 = 0,025 e 15 g.l.; B0_inf = Beta0 - t_tab*s*sqrt(G00); B0_sup = Beta0 + t_tab*s*sqrt(G00); B1_inf = Beta1 - t_tab*s*sqrt(G11); B1_sup = Beta1 + t_tab*s*sqrt(G11); B2_inf = Beta2 - t_tab*s*sqrt(G22); B2_sup = Beta2 + t_tab*s*sqrt(G22); B3_inf = Beta3 - t_tab*s*sqrt(G33); B3_sup = Beta3 + t_tab*s*sqrt(G33); print 'Exemplo 8.6.2',,Beta[format=8.4] s[format=8.4],, 't_tab(0,05/2; 15g.l.)' t_tab[format=8.4] ,,,; print '----------------------------------------------------------------', 'Exemplo 8.6.2: Intervalo de Confiança (gama=95%) para cada Betaj', '----------------------------------------------------------------',; print 'IC(Beta0, 95%) = [' B0_inf[format=8.4] ',' B0_sup[format=8.4] ']'; print 'IC(Beta1, 95%) = [' B1_inf[format=8.4] ',' B1_sup[format=8.4] ']'; print 'IC(Beta2, 95%) = [' B2_inf[format=8.4] ',' B2_sup[format=8.4] ']'; print 'IC(Beta3, 95%) = [' B3_inf[format=8.4] ',' B3_sup[format=8.4] ']'; d = 4; * número de betas; alfa = 0.05/(2*d); t_Bon = tinv(1-alfa,n-k-1); * Calcula t-tabelado para alfa = 0,05/2 = 0,025 e 15 g.l.; B0_inf = Beta0 - t_bon*s*sqrt(G00); B0_sup = Beta0 + t_bon*s*sqrt(G00); B1_inf = Beta1 - t_bon*s*sqrt(G11); B1_sup = Beta1 + t_bon*s*sqrt(G11); B2_inf = Beta2 - t_bon*s*sqrt(G22); B2_sup = Beta2 + t_bon*s*sqrt(G22); B3_inf = Beta3 - t_bon*s*sqrt(G33); B3_sup = Beta3 + t_bon*s*sqrt(G33); print '----------------------------------------------------------------', 'Exemplo 8.6.2: Intervalo de Confiança (gama=95%) para cada Betaj', ' usando o Método de Bonferroni ', '----------------------------------------------------------------',; print 'IC(Beta0, 95%) = [' B0_inf[format=8.4] ',' B0_sup[format=8.4] ']'; print 'IC(Beta1, 95%) = [' B1_inf[format=8.4] ',' B1_sup[format=8.4] ']'; print 'IC(Beta2, 95%) = [' B2_inf[format=8.4] ',' B2_sup[format=8.4] ']'; print 'IC(Beta3, 95%) = [' B3_inf[format=8.4] ',' B3_sup[format=8.4] ']'; quit; * ----------------------- ; * Solução usando proc reg ; * ----------------------- ; proc reg data=ChemReaction; model y1 = x1 x2 x3 / clb; * Opção clb => IC(beta); B1_igual_B2_igual_0: test x1=x2=0; * Testa H0: B1 = B2 = 0; B1_igual_B2: test x1=x2; * Testa H0: B1 = B2; B1_mais_B2_igual_B3: test x1+x2=x3; * Testa H0: B1+B2 = B3; run;