Uma das hipóteses básicas da modelagem linear é constante, ou homogênea, a variância. O que isso significa exatamente? Vamos simular alguns dados que satisfazem esta condição para ilustrar o conceito.
Below criamos um vetor de números ordenado de 1 a 10 chamado x, e depois criamos um vetor de números chamado y que é uma função de x. Quando plotamos x vs y, obtemos uma linha reta com uma intercepção de 1.2 e uma inclinação de 2.1.
x <- seq(1,10, length.out = 100)y <- 1.2 + 2.1 * xplot(x, y)
Agora vamos adicionar algum “ruído” aos nossos dados para que y não seja completamente determinado por x. Podemos fazer isso extraindo aleatoriamente valores de uma distribuição Normal teórica com média 0 e alguma variância definida, e depois adicionando-os à fórmula que gera y. A função rnorm em R permite-nos fazer isto facilmente. Abaixo desenhamos 100 valores aleatórios a partir de uma distribuição Normal com média 0 e desvio padrão 2 e salvamos como um vetor chamado noise. (Lembre-se que o desvio padrão é simplesmente a raiz quadrada da variância.) Depois geramos y com o ruído adicionado. A função set.seed(222) permite obter os mesmos dados “aleatórios” no caso de você querer seguir adiante. Finalmente nós replotamos os dados.
set.seed(222)noise <- rnorm(n = 100, mean = 0, sd = 2)y <- 1.2 + 2.1 * x + noiseplot(x, y)
Agora nós temos dados que é uma combinação de um componente linear, determinístico (\(y = 1,2 + 2,1x\)) e ruído aleatório retirado de uma distribuição \(N(0, 2)\). Estas são as suposições básicas que fazemos sobre os nossos dados quando nos ajustamos a um modelo linear tradicional. Abaixo usamos a função lm para “recuperar” os valores “verdadeiros” usados para gerar os dados, dos quais existem três:
- A intercepção: 1.2
- A inclinação: 2.1
- O desvio padrão: 2
m <- lm(y ~ x)summary(m)## ## Call:## lm(formula = y ~ x)## ## Residuals:## Min 1Q Median 3Q Max ## -5.5831 -1.2165 0.3288 1.3022 4.3714 ## ## Coefficients:## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 1.27426 0.44720 2.849 0.00534 ** ## x 2.09449 0.07338 28.541 < 2e-16 ***## ---## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## ## Residual standard error: 1.926 on 98 degrees of freedom## Multiple R-squared: 0.8926, Adjusted R-squared: 0.8915 ## F-statistic: 814.6 on 1 and 98 DF, p-value: < 2.2e-16A (Intercept) e x (ou seja, a inclinação) as estimativas na secção Coeficientes, 1,27 e 2,09, estão bastante próximas de 1,2 e 2,1, respectivamente. O erro padrão residual de 1,926 também está bem próximo do valor constante de 2. Produzimos um modelo “bom” porque sabíamos como y foi gerado e demos à função lm o modelo “correto” para caber. Embora não possamos fazer isso na vida real, é um exercício útil para nos ajudar a entender as suposições da modelagem linear.
Agora e se a variância não fosse constante? E se multiplicássemos o desvio padrão de 2 pela raiz quadrada de x? Como vemos no gráfico abaixo, a dispersão vertical dos pontos aumenta como x aumentos.
set.seed(222)noise <- rnorm(n = 100, mean = 0, sd = 2 * sqrt(x))y <- 1.2 + 2.1 * x + noiseplot(x, y)
Multiplicamos 2 por sqrt(x) porque estamos especificando o desvio padrão. Se pudéssemos especificar o desvio, teríamos multiplicado 4 por apenas x.
Se cabermos no mesmo modelo usando lm, obtemos um erro Padrão Residual de 4.488.
m2 <- lm(y ~ x)summary(m2)## ## Call:## lm(formula = y ~ x)## ## Residuals:## Min 1Q Median 3Q Max ## -15.0460 -2.4013 0.5638 2.8734 10.2996 ## ## Coefficients:## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 1.313 1.042 1.26 0.211 ## x 2.096 0.171 12.26 <2e-16 ***## ---## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## ## Residual standard error: 4.488 on 98 degrees of freedom## Multiple R-squared: 0.6051, Adjusted R-squared: 0.6011 ## F-statistic: 150.2 on 1 and 98 DF, p-value: < 2.2e-16 Sabemos que isso não está certo, porque simulamos os dados. Não houve desvio padrão constante quando criamos o “ruído”. Cada valor aleatório foi extraído de uma distribuição Normal diferente, cada um com média 0 mas um desvio padrão que variou de acordo com x. Isto significa que a nossa suposição de variância constante é violada. Como detectaríamos isto na vida real?
A forma mais comum é traçar os resíduos versus valores ajustados. Isto é fácil de fazer em R. Basta chamar plot no objeto modelo. Isto gera quatro gráficos diferentes para avaliar as suposições tradicionais de modelagem. Veja este post no blog para mais informações. Os gráficos que nos interessam são o 1º e 3º gráficos, que podemos especificar com o argumento which
plot(m2, which = 1)
plot(m2, which = 3)
No primeiro gráfico vemos a variabilidade em torno de 0 aumentando à medida que nos movemos mais para a direita com valores encaixados maiores. No terceiro gráfico vemos também o aumento da variabilidade à medida que nos movemos para a direita, embora desta vez os resíduos tenham sido padronizados e transformados para a escala da raiz quadrada. A trajetória positiva da linha vermelha lisa indica um aumento na variância.
Então agora que confirmamos que a nossa suposição de variância não-constante é inválida, o que podemos fazer? Uma abordagem é transformar os dados em log. Isto às vezes funciona se a sua variável de resposta for positiva e altamente enviesada. Mas esse não é realmente o caso aqui. y é apenas ligeiramente enviesada. (Chame hist() em y para verificar.) Além disso, sabemos que nossos dados não estão em uma escala de log.
Outra abordagem é modelar a variância não-constante. Esse é o tópico deste blog post.
Para fazer isso, vamos usar funções do pacote nlme que está incluído na instalação base do R. A função workhorse é gls, que significa “generalized least squares” (mínimos quadrados generalizados). Usamo-la tal como usaríamos a função lm, excepto que também usamos o argumento weights juntamente com uma mão cheia de funções de variância. Vamos demonstrar com um exemplo.
Below nós ajustamos o modelo “correto” aos nossos dados que exibiram variância não-constante. Nós carregamos o pacote nlme para que possamos usar a função gls1. Nós especificamos a sintaxe do modelo como antes, y ~ x. Então usamos o argumento weights para especificar a função de variância, neste caso varFixed, também parte do pacote nlme. Isto diz que nossa função de variância não tem parâmetros e uma única covariada, x, que é precisamente como nós geramos os dados. A sintaxe, ~x, é uma fórmula unilateral que pode ser lida como “variância de modelo como uma função de x“
library(nlme)vm1 <- gls(y ~ x, weights = varFixed(~x))summary(vm1)## Generalized least squares fit by REML## Model: y ~ x ## Data: NULL ## AIC BIC logLik## 576.2928 584.0477 -285.1464## ## Variance function:## Structure: fixed weights## Formula: ~x ## ## Coefficients:## Value Std.Error t-value p-value## (Intercept) 1.369583 0.6936599 1.974431 0.0511## x 2.085425 0.1504863 13.857900 0.0000## ## Correlation: ## (Intr)## x -0.838## ## Standardized residuals:## Min Q1 Med Q3 Max ## -2.8942967 -0.6293867 0.1551594 0.6758773 2.2722755 ## ## Residual standard error: 1.925342 ## Degrees of freedom: 100 total; 98 residualO resultado produz um erro Padrão Residual de 1,925 que é muito próximo de 2, o valor “verdadeiro” que usamos para gerar os dados. Os valores (Interceptar) e inclinação, 1,37 e 2,09, também estão muito próximos de 1,2 e 2,1,
Após novamente podemos chamar de plot no objeto modelo. Neste caso apenas um gráfico é gerado: resíduos padronizados versus valores ajustados:
plot(vm1)
Este gráfico parece bom. Desde que modelemos a nossa variância em função de x, o modelo ajustado não se ajusta nem em excesso nem em baixo de forma sistemática (ao contrário de quando usamos lm para ajustar o modelo e assumimos uma variância constante.)
A função varFixed cria pesos para cada observação, simbolizados como \(w_i\). A idéia é que quanto maior o peso de uma determinada observação, menor a variância da distribuição Normal da qual foi extraída a sua componente sonora. No nosso exemplo, como x aumenta a variância. Portanto, pesos mais altos devem ser associados a valores menores de x. Isto pode ser feito tomando o recíproco de x, ou seja \(w_i = 1/x\). Então quando \(x = 2\), \(w = 1/2\). Quando \\(x = 10\), \(w = 1/10\). Maior x obter pesos menores.
Finalmente, para garantir que pesos maiores estão associados a variações menores, dividimos a variação constante pelo peso. Matematicamente,
Or tomando a raiz quadrada para expressar o desvio padrão,
Então quanto maior o denominador (ou seja, quanto maior o peso e, portanto, menor o x), menor a variância e mais preciso o observado y.
Acontece que podemos usar lm para pesar também as observações. Também tem um argumento weights como a função gls. A única diferença é que temos que ser mais explícitos sobre como expressamos os pesos. Neste caso, temos de especificar o recíproco de x. Note que o resultado é quase idêntico ao que obtemos usando gls e varFixed.
m3 <- lm(y ~ x, weights = 1/x)summary(m3)## ## Call:## lm(formula = y ~ x, weights = 1/x)## ## Weighted Residuals:## Min 1Q Median 3Q Max ## -5.5725 -1.2118 0.2987 1.3013 4.3749 ## ## Coefficients:## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 1.3696 0.6937 1.974 0.0511 . ## x 2.0854 0.1505 13.858 <2e-16 ***## ---## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## ## Residual standard error: 1.925 on 98 degrees of freedom## Multiple R-squared: 0.6621, Adjusted R-squared: 0.6587 ## F-statistic: 192 on 1 and 98 DF, p-value: < 2.2e-16 O poder do pacote nlme é que ele permite uma variedade de funções de variância. A função varFixed que acabamos de ilustrar é a mais simples e algo que pode ser feito com lm. As outras funções de variância incluem:
varIdentvarPowervarExpvarConstPowervarComb
Vamos explorar cada uma delas usando dados simulados.
varIdent
A função varIdent permite-nos modelar diferentes variâncias por estrato. Para ver como funciona, vamos primeiro simular os dados com esta propriedade. Abaixo usamos set.seed(11) no caso de alguém querer simular os mesmos dados aleatórios. Depois definimos n igual a 400, o número de observações que vamos simular. x é gerado o mesmo que antes. Neste exemplo, incluímos um preditor adicional chamado g que pode ser pensado como gênero. Nós amostramos aleatoriamente 400 valores de “m” e “f”. A seguir simulamos um vetor de dois desvios padrão e salvamos como msd. Se g == "m", então o desvio padrão é \(2 \ 2,5). Caso contrário, é {\i1}(2}) para g == "f". Usamos este vector na linha seguinte para gerar y. Observe onde msd está plugado na função rnorm. Observe também como nós geramos y. Incluímos uma interação entre x e y. Quando g == "f", a intercepção e a inclinação é 1.2 e 2.1. Quando g == "m", a intercepção e a inclinação são (1.2 + 2.8) e (2.1 + 2.8), respectivamente. Finalmente colocamos nossos dados simulados em um quadro de dados e plotamos com ggplot2.
set.seed(11)n <- 400x <- seq(1,10, length.out = n)g <- sample(c("m","f"), size = n, replace = TRUE)msd <- ifelse(g=="m", 2*2.5, 2)y <- 1.2 + 2.1*x - 1.5*(g == "m") + 2.8*x*(g == "m") + rnorm(n, sd = msd)d <- data.frame(y, x, g)library(ggplot2)ggplot(d, aes(x, y, color = g)) + geom_point()
Notem as diferentes variações para cada nível de g. A variância para “m” é muito maior do que a variância para “f”. Ela tem muito mais dispersão do que “f”. Nós simulamos os dados dessa forma. Definimos a variância para “m” como sendo 2,5 vezes a de “f”.
Agora vamos usar a função gls com varIdent na tentativa de recuperar estes valores verdadeiros. Usamos da mesma forma que antes: definamos nossa função de variância no argumento weights. Abaixo especificamos no argumento form que a fórmula para modelar a variância é condicional a g. A expressão ~ 1|g é uma fórmula unilateral que diz que a variância difere entre os níveis de g. A 1 apenas significa que não temos covariáveis adicionais em nosso modelo para variância. É possível incluir uma fórmula como ~ x|g, mas isso seria incorreto neste caso já que não usamos x na geração da variância. Olhando o gráfico também mostra que enquanto a variabilidade em y é diferente entre os grupos, a variabilidade não aumenta como x aumentos.
vm2 <- gls(y ~ x*g, data = d, weights = varIdent(form = ~ 1|g))summary(vm2)## Generalized least squares fit by REML## Model: y ~ x * g ## Data: d ## AIC BIC logLik## 2081.752 2105.64 -1034.876## ## Variance function:## Structure: Different standard deviations per stratum## Formula: ~1 | g ## Parameter estimates:## f m ## 1.00000 2.58127 ## ## Coefficients:## Value Std.Error t-value p-value## (Intercept) 0.9686349 0.3237724 2.99172 0.0029## x 2.1222707 0.0525024 40.42239 0.0000## gm -1.9765090 0.9352152 -2.11343 0.0352## x:gm 2.9957974 0.1553551 19.28355 0.0000## ## Correlation: ## (Intr) x gm ## x -0.901 ## gm -0.346 0.312 ## x:gm 0.304 -0.338 -0.907## ## Standardized residuals:## Min Q1 Med Q3 Max ## -2.74115039 -0.67013954 0.01619031 0.69793776 2.72985748 ## ## Residual standard error: 2.005397 ## Degrees of freedom: 400 total; 396 residualPrimeiro aviso na parte inferior do output que o erro padrão residual estimado é 2.0053974, muito próximo do valor “verdadeiro” de 2. Também notamos na seção “Função de desvio” que obtemos um valor estimado de 2,58127 para o grupo “m”, que é muito próximo do valor “verdadeiro” de 2,5 que usamos para gerar o desvio diferente para o grupo “m”. Em geral, quando se usa a função varIdent para estimar diferentes variâncias entre níveis de estratificação, um dos níveis será definido como linha de base, e os outros serão estimados como múltiplos do erro padrão residual. Neste caso, como “f” vem antes de “m” em ordem alfabética, “f” foi definido como linha de base, ou 1. O erro padrão residual estimado para o grupo “f” é \ 2.005397 \ 1\). O erro padrão residual estimado para o grupo “m” é \(2.005397 \times 2.58127\).
É importante notar que estes são apenas estimativas. Para ter uma idéia da incerteza dessas estimativas, podemos usar a função intervals do pacote nlme para obter intervalos de confiança de 95%. Para reduzir o output, nós salvamos o resultado e vemos elementos selecionados do objeto.
int <- intervals(vm2) O elemento varStruct contém o intervalo de confiança de 95% para o parâmetro estimado na função de variância. O parâmetro neste caso é o multiplicador de erro padrão residual para o grupo masculino.
int$varStruct## lower est. upper## m 2.245615 2.58127 2.967095## attr(,"label")## "Variance function:" O elemento sigma contém o intervalo de confiança de 95% para o erro padrão residual.
int$sigma## lower est. upper ## 1.818639 2.005397 2.211335 ## attr(,"label")## "Residual standard error:"Os intervalos são bastante pequenos e contêm o valor “verdadeiro” que usámos para gerar os dados.
varPower
A função varPower permite-nos modelar a variância como uma potência do valor absoluto de uma covariada. Mais uma vez, para ajudar a explicar isto, vamos simular os dados com esta propriedade. Abaixo o principal a notar é o argumento sd da função rnorm. É aí que pegamos o desvio padrão de 2 e depois o multiplicamos pelo valor absoluto de x elevado à potência de 1,5. Isto é similar a como simulamos dados para demonstrar a função varFixed acima. Nesse caso simplesmente assumimos que o expoente era 0,5. (A recordação tomando a raiz quadrada de um número é equivalente a elevá-lo para a potência de 0,5). Aqui nós escolhemos arbitrariamente um poder de 1,5. Quando usamos gls com varPower tentaremos recuperar o valor “verdadeiro” de 1,5 além de 2,
set.seed(4)n <- 1000x <- seq(1,10,length.out = n)y <- 1.2 + 2.1*x + rnorm(n, sd = 2*abs(x)^1.5)d <- data.frame(y, x)ggplot(d, aes(x, y)) + geom_point()
Vemos os dados exibindo a forma clássica de variância aumentando à medida que o preditor aumenta. Para modelar corretamente esses dados usando gls, nós os fornecemos com a fórmula y ~x e usamos o argumento weights com varPower. Note que especificamos a fórmula unilateral tal como fizemos com a função varFixed. Neste modelo, no entanto, vamos obter uma estimativa para a potência.
vm3 <- gls(y ~ x, data = d, weights = varPower(form = ~x))summary(vm3)## Generalized least squares fit by REML## Model: y ~ x ## Data: d ## AIC BIC logLik## 8840.188 8859.811 -4416.094## ## Variance function:## Structure: Power of variance covariate## Formula: ~x ## Parameter estimates:## power ## 1.520915 ## ## Coefficients:## Value Std.Error t-value p-value## (Intercept) 2.321502 0.4750149 4.887218 0## x 1.582272 0.2241367 7.059404 0## ## Correlation: ## (Intr)## x -0.845## ## Standardized residuals:## Min Q1 Med Q3 Max ## -2.892319514 -0.655237190 0.001653178 0.691241690 3.346273816 ## ## Residual standard error: 1.872944 ## Degrees of freedom: 1000 total; 998 residualint <- intervals(vm3)int$varStruct## lower est. upper## power 1.445064 1.520915 1.596765## attr(,"label")## "Variance function:"int$sigma## lower est. upper ## 1.650987 1.872944 2.124740 ## attr(,"label")## "Residual standard error:" A potência é estimada em 1,52 o que está muito próximo do valor “verdadeiro” de 1,5. O erro padrão residual estimado de 1,8729439 também está muito próximo de 2. Ambos os intervalos capturam os valores que usamos para simular os dados. Os coeficientes no modelo, por outro lado, são um pouco mal estimados. Isto não é surpreendente dada a variabilidade em y, especialmente para \(x > 2\).
varExp
A função varExp permite-nos modelar a variância como uma função exponencial de uma covariada. Mais uma vez vamos explicar esta função de variância usando dados simulados. A única mudança está no argumento sd da função rnorm. Temos um valor fixo de 2 que multiplicamos por x que por sua vez é multiplicado por 0,5 e exponenciado. Note como a variância aumenta rapidamente à medida que aumentamos x. Isto é devido ao crescimento exponencial da variância.
set.seed(55)n <- 400x <- seq(1, 10, length.out = n)y <- 1.2 + 2.1*x + rnorm(n, sd = 2*exp(0.5*x))d <- data.frame(y, x)ggplot(d, aes(x, y)) + geom_point()
Para trabalharmos ao contrário e recuperarmos estes valores usamos a função varExp no argumento weights de gls. A fórmula unilateral não muda neste caso. Ela diz modelar a variância como uma função de x. A função varExp diz que x foi multiplicada por algum valor e exponenciada, então gls tentará estimar esse valor.
vm4 <- gls(y ~ x, data = d, weights = varExp(form = ~x))summary(vm4)## Generalized least squares fit by REML## Model: y ~ x ## Data: d ## AIC BIC logLik## 3873.153 3889.098 -1932.576## ## Variance function:## Structure: Exponential of variance covariate## Formula: ~x ## Parameter estimates:## expon ## 0.4845623 ## ## Coefficients:## Value Std.Error t-value p-value## (Intercept) 1.198525 1.1172841 1.072713 0.284## x 2.073297 0.4933502 4.202486 0.000## ## Correlation: ## (Intr)## x -0.892## ## Standardized residuals:## Min Q1 Med Q3 Max ## -3.16976627 -0.69525696 0.00614222 0.63718022 2.90423217 ## ## Residual standard error: 2.119192 ## Degrees of freedom: 400 total; 398 residualint <- intervals(vm4)int$varStruct## lower est. upper## expon 0.4562871 0.4845623 0.5128374## attr(,"label")## "Variance function:"int$sigma## lower est. upper ## 1.786737 2.119192 2.513505 ## attr(,"label")## "Residual standard error:" A estimativa “expon” de 0,4845623 na seção “função de variância” é muito próxima do nosso valor especificado de 0,5. Da mesma forma, o erro padrão residual estimado de 2,1191918 está próximo do valor “verdadeiro” de 2. As estimativas do coeficiente do modelo também estão próximas dos valores que usamos para gerar os dados, mas observe a incerteza no (Interceptar). O teste de hipóteses no resumo não pode descartar uma intercepção negativa. Isto novamente não é surpreendente já que y tem tanta variância não-constante assim como o fato de não termos observações de x igual a 0. Como a intercepção é o valor y toma quando x igual a 0, nossa intercepção estimada é uma extrapolação para um evento que não observamos.
varConstPower
A função varConstPower permite-nos modelar a variância como uma constante positiva mais uma potência do valor absoluto de uma covariada. Isso é uma boca cheia, mas isto é basicamente o mesmo que a função varPower, exceto que agora adicionamos uma constante positiva a ela. O seguinte código simula dados para os quais a função varConstPower seria adequada para ser usada. Note que é idêntico ao código que usamos para simular dados para a seção varPower, exceto que adicionamos 0,7 a x na função abs. Porquê 0.7? Nenhuma razão especial. É apenas uma constante positiva que escolhemos.
set.seed(4)n <- 1000x <- seq(1,10,length.out = n)y <- 1.2 + 2.1*x + rnorm(n, sd = 2*abs(0.7 + x)^1.5)d <- data.frame(y, x)ggplot(d, aes(x, y)) + geom_point()
> A forma correcta de modelar estes dados é usar gls com a função varConstPower. A fórmula unilateral é a mesma de antes. O resumo retorna três estimativas para o modelo de variância: a constante, a potência e o erro padrão residual. Os valores “verdadeiros” são 0,7, 1,5 e 2, respectivamente.
vm5 <- gls(y ~ x, data = d, weights = varConstPower(form = ~x))summary(vm5)## Generalized least squares fit by REML## Model: y ~ x ## Data: d ## AIC BIC logLik## 9319.794 9344.323 -4654.897## ## Variance function:## Structure: Constant plus power of variance covariate## Formula: ~x ## Parameter estimates:## const power ## 0.3554921 1.3593016 ## ## Coefficients:## Value Std.Error t-value p-value## (Intercept) 2.763788 0.7607449 3.633003 3e-04## x 1.456677 0.2915621 4.996112 0e+00## ## Correlation: ## (Intr)## x -0.824## ## Standardized residuals:## Min Q1 Med Q3 Max ## -2.846999539 -0.664551194 0.006709051 0.679895874 3.287366683 ## ## Residual standard error: 2.887761 ## Degrees of freedom: 1000 total; 998 residualint <- intervals(vm5)int$varStruct## lower est. upper## const 0.02658804 0.3554921 4.753063## power 1.14258677 1.3593016 1.576016## attr(,"label")## "Variance function:"int$sigma## lower est. upper ## 1.814757 2.887761 4.595195 ## attr(,"label")## "Residual standard error:"Os intervalos capturam os valores “verdadeiros”, embora o intervalo para a constante seja bastante amplo. Se não soubéssemos o valor verdadeiro, pareceria que a constante poderia plausivelmente ser 2 ou mesmo 4,
varComb
Finalmente, a função varComb permite-nos modelar combinações de dois ou mais modelos de variância multiplicando em conjunto as funções de variância correspondentes. Obviamente isto pode acomodar alguns modelos de variância muito complexos. Vamos simular alguns dados básicos que seriam apropriados para usar com a função varComb.
O que fazemos abaixo é combinar dois processos de variância diferentes:
- Um que permite que o desvio padrão seja diferente entre os sexos (“f” = \(2\), “m” = \(2 \ vezes 2.5\))
- Um que permite que o desvio padrão aumente como
xaumenta, ondexé multiplicado por 1,5 e exponenciado
Para ajudar a visualizar os dados que temos limitado x a variar de 1 a 3.
set.seed(77)n <- 400x <- seq(1, 3, length.out = n)g <- sample(c("m","f"), size = n, replace = TRUE)msd <- ifelse(g=="m", 2*2.5, 2) * exp(1.5*x)y <- 1.2 + 2.1*x - 1.5*(g == "m") + 2.8*x*(g == "m") + rnorm(n, sd = msd)d <- data.frame(y, x, g)ggplot(d, aes(x, y, color = g)) + geom_point() 
O gráfico mostra um aumento da variância como x aumenta, mas também as diferenças de variância entre os sexos. A variância de y para o grupo “m” é muito maior que a variância de y no grupo “f”, especialmente quando x é maior que 1,5,
Para modelar corretamente o processo de geração de dados que especificamos acima e tentar recuperar os valores verdadeiros, usamos a função varComb como um invólucro em torno de mais duas funções de variância: varIdent e varExp. Por que estas duas? Porque temos variâncias diferentes entre gêneros, e temos variância aumentando exponencialmente como uma função de x.
vm6 <- gls(y ~ x*g, data = d, weights = varComb(varIdent(form = ~ 1|g), varExp(form = ~x)))summary(vm6)## Generalized least squares fit by REML## Model: y ~ x * g ## Data: d ## AIC BIC logLik## 4431.45 4459.32 -2208.725## ## Combination of variance functions: ## Structure: Different standard deviations per stratum## Formula: ~1 | g ## Parameter estimates:## f m ## 1.000000 2.437046 ## Structure: Exponential of variance covariate## Formula: ~x ## Parameter estimates:## expon ## 1.540608 ## ## Coefficients:## Value Std.Error t-value p-value## (Intercept) -5.996337 6.539552 -0.9169339 0.3597## x 8.829971 4.849272 1.8208860 0.0694## gm -17.572238 16.548540 -1.0618603 0.2889## x:gm 16.932938 12.192793 1.3887661 0.1657## ## Correlation: ## (Intr) x gm ## x -0.972 ## gm -0.395 0.384 ## x:gm 0.387 -0.398 -0.974## ## Standardized residuals:## Min Q1 Med Q3 Max ## -2.74441479 -0.67478610 -0.00262221 0.66079254 3.14975288 ## ## Residual standard error: 1.794185 ## Degrees of freedom: 400 total; 396 residualA seção de resumo contém duas seções para modelagem de variância. A primeira estima o multiplicador para o grupo “m” em 2,437, que está muito próximo do valor real de 2,5. O parâmetro exponencial é estimado em 1,54, extremamente próximo do valor verdadeiro de 1,5 que usamos na geração dos dados. Finalmente o erro padrão residual é estimado em cerca de 1,79, próximo ao valor verdadeiro de 2. A chamada intervals(vm6) mostra intervalos de confiança muito apertados. Os prefixos A e B no elemento varStruct são apenas etiquetas para os dois modelos de variância diferentes.
int <- intervals(vm6)int$varStruct## lower est. upper## A.m 2.119874 2.437046 2.801673## B.expon 1.419366 1.540608 1.661850## attr(,"label")## "Variance function:"int$sigma## lower est. upper ## 1.380008 1.794185 2.332669 ## attr(,"label")## "Residual standard error:"Felizmente devido à grande variabilidade exponencial, as estimativas dos coeficientes do modelo são terrivelmente ruins.
Para quê?
Então porque é que simulámos todos estes dados falsos e depois tentamos recuperar valores “verdadeiros”? De que serve isso? Perguntas justas. A resposta é que nos ajuda a compreender que suposições estamos a fazer quando especificamos e encaixamos um modelo. Quando especificamos e encaixamos o seguinte modelo…
m <- lm(y ~ x1 + x2)…estamos a dizer que pensamos que y é aproximadamente igual a uma soma ponderada de x1 e x2 (mais uma intercepção) com os erros a serem extraídos aleatoriamente de uma distribuição Normal com média de 0 e algum desvio padrão fixo. Da mesma forma, se especificarmos e ajustarmos o seguinte modelo…
m <- gls(y ~ x1 + x2, data = d, weights = varPower(form = ~x1))…estamos dizendo que pensamos que y é aproximadamente igual a uma soma ponderada de x1 e x2 (mais uma intercepção) com erros sendo sorteios aleatórios de uma distribuição Normal com média 0 e um desvio padrão que cresce como um múltiplo de x1 levantado por alguma potência.
Se pudermos simular dados adequados para esses modelos, então temos uma melhor compreensão das suposições que estamos fazendo quando usamos esses modelos. Esperemos que agora você tenha um melhor entendimento do que você pode fazer para variar os modelos usando o pacote nlme
Para perguntas ou esclarecimentos sobre este artigo, contate a biblioteca UVA StatLab: [email protected]
Veja toda a coleção de artigos da biblioteca UVA StatLab.
Clay Ford
Statistical Research Consultant
University of Virginia Library
April 7, 2020
-
1. O pacote
nlmeé talvez mais conhecido pela sua funçãolme, que é usada para encaixar modelos de efeito misto (ou seja, modelos com efeitos tanto fixos como aleatórios). Este post do blog demonstra funções de variância usandogls, que não se encaixa em efeitos aleatórios. Entretanto, tudo que apresentamos neste post do blog pode ser usado comlme.