title: Is my parameter estimate significant or not?

text: DISCLAIMER: This question was sent to Stata list today, but so far nobody has answered. NOTE: I use Stata here, but actually I don't think the question is software- specific. Hi, I would appreciate someone could provide me an answer to the following question: I am estimating the following model:               . areg beta L.lev group#cL.lev i.year, absorb(group)      Group is a categorical variable: 1, 2 Lev is a continuous variable Year are year effects Group are group effects Under that specification of the model, I get the following results:                   . areg beta L.lev group#cL.lev i.year, absorb(group)          Linear regression, absorbing indicators           Number of obs   =        285                                                       F(  17,    266) =       4.04                                                       Prob > F        =     0.0000                                                       R-squared       =     0.3540                                                       Adj R-squared   =     0.3103                                                       Root MSE        =     0.3038          ------------------------------------------------------------------------------             beta |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]     -------------+----------------------------------------------------------------              lev |              L1. |    .059685    .014907     4.00   0.000     .0303342    .0890358                  |     group#cL.lev |             2  |  -.0451045   .0178834    -2.52   0.012   -.0803155   -.0098936      What I would like to highlight from here, is the fact that the P value for group#cLlev 2 is significant at the 5% level. If I run the same model in a different way (not including the variable Lev by itself), like this:               . areg beta group#cL.lev i.year, absorb(group)          Linear regression, absorbing indicators           Number of obs   =        285                                                       F(  17,    266) =       4.04                                                       Prob > F        =     0.0000                                                       R-squared       =     0.3540                                                       Adj R-squared   =     0.3103                                                       Root MSE        =     0.3038          ------------------------------------------------------------------------------             beta |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]     -------------+----------------------------------------------------------------     group#cL.lev |               1  |    .059685    .014907     4.00   0.000     .0303342    .0890358               2  |   .0145805   .0114004     1.28   0.202    -.0078661     .037027      The P value for group 2 is not significant at the 5% level. So the question is, which P-value should I consider as correct? In other words, is my parameter estimate significant or not? (Please note that I get exactly the same parameter estimates from both methods, but the associated P-values and t-values change).


title: Rolling out trial in phases, but not stepped wedge

text: **Update (Sept 10, 2013)** : I believe it would be more correct to say that increasing baseline or endline measurements are ways to decrease the design effect, thereby making the SW design efficient, rather than stating that the maximum number of measurements is necessary. Woertman et al. (2013). The stepped wedge design (pdf) is a nice alternative to parallel-group designs when, for logistical reasons, the intervention must be rolled-out in stages. A potential downside of this design, however, can be the number of measurement rounds. Even though SW designs can have increased power (thus reducing required sample size to detect the same effect), every unit is observed/measured before and after every treatment round (step). If you have five steps, there are six measurement rounds, including the initial measurement round when all units are in the control group. So if you have n=1000, this is 1000 x 6 = 6000 observations/measurements. I'm writing to ask about a possible alternative (see pic below):   1. Stratifying by community (let's say 5 communities overall; increase to N = 1500 because design has less power than SW)   2. Randomize strata (communities) to order of intervention (first, second, third, fourth, fifth)   3. Within the first community strata, Community A, conduct baseline surveys with all n=300 and then randomize units to treatment or control   4. Deliver the intervention to n/2 units randomized to treatment   5. Conduct endline survey with all n=300 in Community A (treatment and wait-list control)   6. Conduct baseline survey with all n=300 in Community B (could be at same time as #5) and then randomize units to treatment or control   7. Deliver the intervention to n/2 units in Community B randomized to treatment AND n/2 wait-list control units from Community A (optional, but this is what we would do)   8. Repeat. In the alternative design, every unit is surveyed twice, just at different times. With a total sample of n=1500, this is 1500 x 2 = 3000 surveys. Compared to the SW design, this is 6000 - 3000 = 3000 fewer surveys, which has big cost implications. SW works because we observe every unit before and after every step and then model time. In the alternative design, we only have 2 measurements (baseline and endline) for every unit assigned to treatment (n=750) and wait-list control (n=750). In alternative:   * Baseline for Community A conducted in month 1    * Endline for Community A conducted in month 3    * Baseline for Community B conducted in month 3   * Endline for Community B conducted in month 6    * Baseline for Community C conducted in month 6    * Endline for Community C conducted in month 9   * Baseline for Community D conducted in month 9    * Endline for Community D conducted in month 12    * Baseline for Community E conducted in month 12   * Endline for Community E conducted in month 15   * (would not measure post-treatment for Community E wait-list control; just deliver program) In the alternative design, can we account for the fact that the observations are made at different times? In SW, every unit is measured before and after every round, which makes it easier to model time effects. Could we regress endline DVs on assignment to treatment (0/1), a vector of baseline controls, dummy variables for community strata, and month of endline measurement? Better alternatives? Assuming there is a solution, how to think through the implications for power? Alternative design: ![enter image description here](http://i.stack.imgur.com/AFdrZ.jpg)


title: Mediation and moderation in labeling

text: So I'm looking at student stigma affecting self-esteem. The overall idea is that stigma affects student achievement through the mediator self-esteem, so I will do multiple regressions and test for mediation. The issue here is that I want to look at moderators. I know that moderators can be gender and ethnicity. My question is, how do I test for moderation in this case? Do regressions with stigma and achievement by gender and ethnicity?


title: Working with residuals of regression

text: So the background is that the I collected yield data for past 5-6 decades and location from where I collected yield data had high yielding varieties introduced over time. I am looking at the relationship between yield and rainfall but this introduction of HYV might affect the true impact of monsoon on yield and therefore I am detrending the data to remove the effect of HYV. I did a linear regression of yield against time in R:               mdl1 <- lm(yield ~ time, data=data)      and then removed the linear trend by taking the residuals of the above regression:               yield.res <- resid(mdl1)      Now I am using these residuals for my subsequent analysis. For example, the relationship between yield and rainfall is:                mdl2 <- lm(yield.res ~ rain, data=data)      In this case, do my `yield.res` have to be normally distributed before I do this regression? If yes, what sort of transformation do I need to use? Since `yield.res` consists of both negative and positive numbers, I am slightly confused how to go about it.


title: Time-series stationarity

text: If I difference a time series and take out trend and seasonality ... does it mean we are left with only irregularity on which we plot the acf and pacf to arrive at the MA and AR order? Do 1st difference, 2nd difference always detrend the series, or do we need to detrend separately?


title: Estimating probabilities of probability products

text: I have a process which generates a binary outcome from a multiplication of two independent stochastic binary inputs: Out = X1(P1) * X2(P2) . X1 and X2 generate a binary output with some probability, P1 and P2, respectively. I have a set of examples, for which I know "Out" and P1. I would like to estimate P2. Note: P2 is constant across all examples, but P1 and Out are not. What is a good way to do this? Thanks!


title: why minimize loss function instead of maximizing reward function?

text: Why is the "de-facto" in statistics to minimize the sum of squared errors cost function instead of maximizing some reward function like the likelihood function?


title: Logarithmic probabilty distribtion and respective problems in matlab

text: Dear crossvalidated community, quite frankly, I'm not that much into statistics and therefore, I'm having much trouble solving some issues im having in Matlab right now. My plan is to imitate/simulate some realistic driving patterns of costumer in a car sharing simulator. Therefore, I picked a logarithmic probability distribution that represents/models the desire to rent a car once in a while. So as of now, I have 500 customer that are devided into three different types of customers with different mean intervals:               if  hourOfDay > 20 || hourOfDay < 9 % night                          if obj.custType == 1                         meanInterval = 45*24*60;                         stDev = meanInterval/10;                         MU = log(meanInterval^2 / sqrt(stDev^2+meanInterval^2));                         SIGMA = sqrt(log(stDev^2/meanInterval^2 + 1));                         obj.minInterval = lognrnd(MU,SIGMA);                          elseif obj.custType == 2                         meanInterval = 30*24*60;                         stDev = meanInterval/8;                         MU = log(meanInterval^2 / sqrt(stDev^2+meanInterval^2));                         SIGMA = sqrt(log(stDev^2/meanInterval^2 + 1));                         obj.minInterval = lognrnd(MU,SIGMA);                          elseif obj.custType == 3                         meanInterval = 25*24*60;                         stDev = meanInterval/9;                         MU = log(meanInterval^2 / sqrt(stDev^2+meanInterval^2));                         SIGMA = sqrt(log(stDev^2/meanInterval^2 + 1));                         obj.minInterval = lognrnd(MU,SIGMA);                     end                      else %daytime                     if obj.custType == 1                         meanInterval = 7*24*60;                         stDev = meanInterval/5;                         MU = log(meanInterval^2 / sqrt(stDev^2+meanInterval^2));                         SIGMA = sqrt(log(stDev^2/meanInterval^2 + 1));                         obj.minInterval = lognrnd(MU,SIGMA);                          elseif obj.custType == 2                         meanInterval = 20*24*60;                         stDev = meanInterval/9;                         MU = log(meanInterval^2 / sqrt(stDev^2+meanInterval^2));                         SIGMA = sqrt(log(stDev^2/meanInterval^2 + 1));                         obj.minInterval = lognrnd(MU,SIGMA);                          elseif obj.custType == 3                         meanInterval = 15*24*60;                         stDev = meanInterval/9;                         MU = log(meanInterval^2 / sqrt(stDev^2+meanInterval^2));                         SIGMA = sqrt(log(stDev^2/meanInterval^2 + 1));                         obj.minInterval = lognrnd(MU,SIGMA);                     end      the "obj.mininterval" is the key factor here (determined by the logarithmic prob. density function) and is used to determine whether they want to drive or not:               hasNotDriven = customers{cc}.minInterval < tt - customers{cc}.lastNightDrive;      (in case that the mininterval is lower than the current time of the simulation minus their last rental, they start driving) For example in this code, customer type 1 is supposed to drive once a week (plus/minus a little deviation so that it's a more or less random pattern) at daytime and, if it is nighttime, each 45 days. In order to avoid that all customers start driving immediately after I started the simulation, I randomly defined when they drove the last time (e.g.):               if obj.custType == 1                     obj.lastDayDrive = -(4*24*60)*rand(1)*rand(1);                     obj.lastNightDrive = -(42.234*24*60)*rand(1)*rand(1);      However, they don't abide by my rules and keep driving and renting cars like it's their sole purpose in life (they drive far more often than just once a week). Have a look at the following image: ![http://imgur.com/6H5UyqM](http://i.stack.imgur.com/w6Jxj.jpg) the green parts of the figure represent daytime and the black ones nighttime. I'm well aware that my request is rather confusing but please feel free to ask as much as you require. I mean I clearly defined the mean, so why do I keep having these spikes?! I mean, as you can see it's more like all 500 customers drive once a week although, from the various mean intervals, they should be driving far less often (solely customer type 1 is supposed to drive once a week, the other customer types should drive like once every twenty days or so; thus, I consider the mininterval to be lower than it actually should be). edit: lowering the standard deviation didn't help either: ![enter image description here](http://i.stack.imgur.com/XK4lC.jpg)


title: Interaction between continuous variable and categorical variable (2+ categories)

text: I'm about to do my analysis of my Master's thesis, but I have a question. I want to check whether media effects on political trust are different between natives and immigrants. These are my variables:   * Dependent: political trust (a continuous variable, ranging from 0 to 70 (I summed up seven items on an 11-point scale, after doing factor analysis)).   * Independent variables:      * TV news: the amount of news people watch on television (11-point scale ranging from 0 (never) to 10 (more than 3,5 hours)).     * TV entertainment: the same for entertainment content.     * Newspaper news     * Newspaper entertainment     * Radio news     * Radio entertainment I created five dummy variables for migrant status: a respondent can either belong to the native group, be a first generation Western immigrant, second generation western immigrant, first generation non-Western immigrant and a second generation non-Western immigrant. This means I have to include four dummies in the analysis. Now I want to check whether the effect of, for example, television news on political trust is different between natives or immigrants. How do I have to make interactions between this continuous variable (TV news consumption) and the categorical variable (migrant status, the five dummies). How many interactions do I have to make? Vnews*the 5 migrant status dummies?


title: How to check the distribution of the given data

text: I am working on time series. I have a set of data which I would like to use for estimation. Can some one tell me how to find under what distribution the data I have goes in. I tried plotting using               histfit(data) in matlab      can some one help. Thanks


title: Pearson correlation for 'lumped' populations?

text: I'm currently using the Pearson (product-moment) correlation coefficient (PPCC) to look at the relationship between deprivation and healthcare usage in a population, but I think I've probably got this wrong. The problem is that the population (about a million people) is broken up into about 100 geographical groups, with group sizes of between 4,000 and 20,000 people. Each group has a single deprivation ($\mathit{x}$) value, and the healthcare usage of the people in the group is combined into a single $\mathit{y}$ value for that group. What I'm doing is taking my 100 ($\mathit{x,y}$) points, and calculating the PPCC for those points. I think (or thought) that this makes sense, because I'm only interested in comparing the geographical groups against each other; I'm not really interested in the individual people per se. The problem with this approach is that a group of 4,000 people is effectively given the same weight as a group 20,000 people. I'm not yet convinced that this is a problem but, if so, is there a way to handle this? Note that I can't just calculate $\rho$ for the ungrouped population of 1,000,000 people, because a single person doesn't have a quantifiable $\mathit{y}$ value. 'Healthcare usage' only makes sense in the context of a group (a group may have 10 healthcare events per 10,000 people per month, for example). One approach that I'm thinking about is simply to exclude outliers. This sort- of makes sense because my detection of outliers takes into account the size of the group to find out whether or not the healthcare usage of the group is within 2$\sigma$ of the entire population mean. Smaller groups are more likely to be outliers, and the $\mathit{(x,y)}$ points that are left after this process can be considered to be valid, irrespective of group size.


title: how to find the probability

text: I transmitted 50 packets to receiver but few of them are missing and i assumed that missing probability of each packet is 0.1.Now i retransmit packet like- 1 0 1 0 1 1 1 0 1 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ------1st retarnsmission row 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 1 0 1 1 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0-------2nd retransmission row 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 1 1 1 1 0 -------3rd retransmission row 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 0 1 1 -------4th retransmission row 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 1 0 1 1 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0-------5th retranmsission row 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 1 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1-------6th retransmission row *1 means packet and 0 means no packet For every row i transmit 9 packets and this is fixed and for 1st 5 row there is no duplicate packet[it means for retransmission purpose this packets are 1st time transmitted] but in six row first 5 packet is original but last 4 packets are already transmitted and this is the second time as i am considering 50 transmitted packet.for 7th row-all 9 packets are duplicate of previous. For 1st 5 retransmission i found that receiver demodulated those packets with probability =0.3887[using bernouli equation-9C1*0.1*(1-0.1)^9-1].6th row demodulation probability is =0.369 as 4 of them are duplicate and 7th row demodualation percentage =0.29 as 9 of them are duplicate.So every case probability is decreasing.but i dont know how this demodulation value comes But i only know in 6th row the missing probability of 4 packet is decreased as it is duplicate and this probability is around 0.049.for 6th row the thing is like- (0.1*5+0.0494*4)/9=0.0775 and demodulation=(9C1*0.0775*(1-0.0775)^9-1=0.365). Now my problem is i cannot find the missing probability of 4 packets of row 6 and also the missing probability of 9 packets of row 7 which is less than 0.1 .How can i calculate that probability?


title: R statistics output interpretation : ebook

text: I am basically a biologist working in genomic data analysis. I have started learning statistics and now am acquainted with basics of it. I can run several statistics commands on R. I have plenty of pdfs which emphasise how to run R or how to perform statistical tests etc., but do not explain the output in a comprehensive way. I am looking for a pdf for R statistics output interpretation, where the output is interpreted in a clear and lucid way. Can anybody share an ebook or pdf having the interpretation of various statistical tests / models?


title: GPML gives too large length scale when optimising hyperparameters

text: I recently started trying to apply Gaussian process regression to a problem, using the MATLAB GPML toolbox. The problem has five (or more) input variables, but for now I'm just looking at one of them. With manually selected hyperparameters, I get reasonable-looking output, but my problem is that when I `minimize()` the hyperparameters, I get a very large length scale, thus producing a flat mean and covariance function. If I run the following code, `hyp2.cov = [5.4, 1.1]` after the call to minimize, so the length scale is `exp(5.4)=214`, even though the data points are on the interval [-0.14, 0.06]. Could anyone tell me why this is so?               x=[-0.0350; 0.0550;-0.0830;-0.1360; 0.0190];     y=[3.0700;3.3200;3.0400;3.0600;2.9200];     covfunc = @covSEiso; hyp2.cov = [0; 0]; hyp2.lik = log(0.1);     likfunc = @likGauss; sn = 0.1; hyp2.lik = log(sn);     hyp2 = minimize(hyp2, @gp, -100, @infExact, [], covfunc, likfunc, x, y);


title: Lasso-ing the order of a lag?

text: Suppose I have longitudinal data of the form $\mathbf Y = (Y_1, \ldots, Y_J) \sim \mathcal N(\mu, \Sigma)$ (I have multiple observations, this is just the form of a single one). I'm interested in restrictions on $\Sigma$. An unrestricted $\Sigma$ is equivalent to taking $$ Y_j = \alpha_j + \sum_{\ell = 1} ^ {j - 1} \phi_{\ell j} Y_{j-\ell} + \varepsilon_j $$ with $\varepsilon_j \sim N(0, \sigma_j)$. This is typically not done since it requires estimating $O(J^2)$ covariance parameters. A model is "lag-$k$" if we take $$ Y_j = \alpha_j + \sum_{\ell = 1} ^ k \phi_{\ell j} Y_{j - \ell} + \varepsilon_j, $$ i.e. we only use the preceding $k$ terms to predict $Y_j$ from the history. What I'd really like to do is use some kind of shrinkage idea to zero out some of the $\phi_{\ell j}$, like the LASSO. But the thing is, I also would like the method I use to prefer models which are lag-$k$ for some $k$; I'd like to penalize higher order lags more than lower order lags. I think this is something we would particularly like to do given that the predictors are highly correlated. An additional issue is that if (say) $\phi_{35}$ is shrunk to $0$ I would also like it if $\phi_{36}$ is shrunk to $0$, i.e. the same lag is used in all of the conditional distributions. I could speculate on this, but I don't want to reinvent the wheel. Is there any LASSO techniques designed to get at this sort of problem? Am I better off just doing something else entirely, like stepwise inclusion of lag orders? Since my model space is small, I could even use an $L_0$ penalty on this problem I guess?


title: Entropy and Likelihood Relationship

text: This is a theoretical question. Suppose that I have a sample s1 coming from distribution K, and a sample s2 coming from distribution M. But I don't know what K or M are. I hypothesize that s1 and s2 are coming from distribution T. I plug in s1 and s2 to pdf of T and calculate their likelihoods, say L(s1), L(s2). If I know differential entropies of K and M, being H(K) and H(M) respectively, can I find a relation between L(s1) and L(s2) in terms of H(K) and H(M). Basically what I am trying to show is something along the lines of this: $\frac {L(s1)}{L(s2)} \sim \frac {e^{H(K)}} {e^{H(M)}}$ Is this even meaningful? Thanks


title: Sampling Technique?

text: I am writing a report in which I have taken the incidents of bomb-blasts as the population and as the sample, I have taken the incidents in which there were some fatalities. kindly tell me which type of non-probability sampling is this? Is it Judgmental sampling or Purposive sampling in which the researcher chooses the sample based on who they think would be appropriate for the study as this is used primarily when there is a limited number of people that have expertise in the area being researched. Also, Kindly tell me what would be my defense if I am asked to justify my sample size? The ume-sakran's table is of no use here. Kindly reply asap. Thank You.


title: Boltzmann machines and the asymmetry between 0 and 1

text: I spent some time yesterday learning about Boltzmann machines (and their cousins, restricted Boltzmann machines). They seem interesting and very powerful, but one little thing is bugging me. The energy function for a Boltzmann machine is (assuming the use of a bias node) $$ E = -\sum_{ij} w_{ij} s_i s_j, $$ where the $w_{ij}$ are the weights, and the $s_i$ are the states of the nodes. Every resource I found used $0$ and $1$ as the possible values for $s_i$, but this seems to make the energy function asymmetric, in the sense that two nodes $i$ and $j$ _both_ have to be in the $1$ state for the $w_{ij}$ term to enter into the energy expression, whereas if just one of them is on, it's the same as them both being off. But it seems that there's nothing in the derivation of the learning algorithm that requires this. As far as I can see, one could use $-1$ and $1$ instead and the learning algorithm would be exactly the same. In this case the contribution to the energy would depend only on whether the two nodes are in the same state or different ones. My question is whether this makes any difference, or whether the two possibilities turn out to be equivalent somehow. If they are different, is it just that using $0$ and $1$ turns out to be better in practice? I can imagine this being the case if Boltzmann machines are typically used in cases where $0$ and $1$ have semantically different meanings, but I find it surprising that this issue isn't mentioned in any of the introductory articles I found online.


title: Distribution of the sum of two independent Beta-Binomial variables

text: Consider two independent discrete random variables $y_1$ and $y_2$, both distributed with a Beta-Binomial distribution, with different number of successes $n_1$ and $n_2$ but the same parameters $a$ and $b$ $ p(y_1|n_1,a,b) = {n_1 \choose y_1} \dfrac{B(y_1 + a,n_1 -y_1 +b)}{B(a,b)} $ $ p(y_2|n_2,a,b) = {n_2 \choose y_2} \dfrac{B(y_2 + a,n_2 -y_2 +b)}{B(a,b)} $ Consider a discrete variable $Z = y_1 + y_2$. Is $Z$ distributed as well as a Beta-Binomial (with parameters $n_1+n_2$, $a'$ and $b'$? I could not prove it in an analytic form so far, but I have been trying it out with some simulations, at least to check whether the assumption is wrong in some cases. Reparametrising $a$ and $b$ as $\mu = \dfrac{a}{a+b}$ and $\rho = \dfrac{1}{a+b+1}$, $y_1$ and $y_2$ have mean $\mu n_1$ and $\mu n_2$ respectively and variance $\mu(1-\mu)n_1(1 + (n_1-1)\rho)$ and $\mu(1-\mu)n_2(1 + (n_2-1)\rho)$ respectively. Based on independence, the mean of $Z$ is $\mu(n_1+n_2)$ and the variance of $Z$ is $\mu(1-\mu)(n_1(1 + (n_1-1)\rho) + n_2(1 + (n_2-1)\rho))$. If $Z$ was distributed according to a beta-binomial distribution, then it would have paramters $\mu'$ and $\rho'$, with $\mu'= \mu$ and $\rho' = \dfrac{\dfrac{n_1(1 + (n_1-1)\rho) + n_2(1 + (n_2-1)\rho)}{n_1+n_2} - 1}{n_1+n_2-1} = \rho \dfrac{n_1(n_1-1)+n_2(n_2-1)}{(n_1+n_2)(n_1+n_2 -1)} $ Here is some code to generate $Z$ as a sum of two independent Beta-Binomials (sorry about the code, R is not my main language)               n1 = 20     n2 = 50     mu = .6     k  = 20          p1  = rbeta(1e6,mu*k,(1-mu)*k)     y1  = rbinom(1e6,n1,p1)          p2  = rbeta(1e6,mu*k,(1-mu)*k)     y2  = rbinom(1e6,n2,p2)          z   = y1+y2          rho  = 1/(k+1)     rho1 = rho*(n1*(n1-1)+n2*(n2-1))/((n1+n2)*(n1+n2-1))     k1   = 1/ rho1 - 1     p3  = rbeta(1e6,mu*k1,(1-mu)*k1)     z1  = rbinom(1e6,n1+n2,p3)          print(c(var(z),var(z1)))     plot(density(z,width= 3))     lines(density(z1,width = 3))      I have been trying this code for different values of $n_1$, $n_2$, $\mu$ and $k$, but in all the cases the variances using the sum of two beta-binomials or an appropriately tuned beta-binomial are very similar (the densities look indistinguishable)


title: Can I Use a Dummy Variable in Quadratic form in lm() in R?

text: I'm trying to replicate models published in The Age-Productivity Gradient: Evidence from a Sample of F1 Drivers and am stumbling at first hurdle over a line "the baseline specification contains a quadratic in age and drivers dummies". I can get quadratic `age` (a numerical var) and linear `name` dummy (a factor in the original dataset) with the following construction:               lm(pos ~ name + I(age) +I((age)^2),data=ergastData)      but I get an error on               lm(pos ~ name + I(name^2) + I(age) +I((age)^2),data=ergastData)      because `name` is a factor and squaring it is not meaningful? So how can I get a quadratic term for the dummy variable? _(I'm not a statistician.. I assume this is a meaningful thing to do and my interpretation of the original paper is correct!)_ As an extension of this question, with ggplot2, I can chart the individual `name` model results using:               g = ggplot(egastdata, aes(x=age,y=pos,col=name))     g = g + stat_smooth(method = "lm", formula = y ~ I(x) +I((x)^2)) +geom_point()     g      This seems to group on `name` and presumably generate separate models for each driver (i.e. a set of distinct `age` coefficients for each driver). How would I use `ggplot` to chart the result of a single set of `age` coefficients and then separate `name` dummy variable coefficients (both linear in `name` and quadratic)? (That is, I want a single model around `age` with a separate coefficients around `name` for each driver.) Would the age coefficients be the same as if I ran the model over the whole set of driver data without distinguishing drivers?


title: Is it correct to say that p-values for Pearson's r decrease with increasing sample size because bias decreases with increased sample size?

text: I realise there are other reasons why p-values will decrease with increased sample size, but I am wondering if this reason is valid at all.


title: Using ordinal regression to evaluate predictor "importance"?

text: We've got a construct-likert-scale with an internal (8 items) and an external dimension (6 item) and there is also a 5-point item y assessing the "subjective" perception ( _How skilled do you think you are?_ ). $N$ is 250 y-Cell frequency:               extrem gering   sehr gering        gering        mittel          hoch     sehr hoch  extrem hoch      3               2                  2             36              78       103            26       Now we've got the hypothesis/idea that subjective perception of the construct is "more" related to the external dimension that to the internal one. In order to prove this, we used three proportional odds regressions (a: `y ~ int + ext`, b: `y ~ int`, c: `y ~ ext`) using the R-package `rms`. Result:   * Every LR-ChiSq p-Value is significant    * Every Intercept and every predictor is highly significant   * The Nakerkes' pseudo-$R^2$ differs      * .23 for a)     * .11 for b)     * .22 for c)  So we used the different $R^2$-values to say "using ext "explains" more than int and including int does only minimal effects once ext is included". What do you think of this approach? Is this even correct?


title: Calculation of an "unconstrained" normal distribution (starting from a censored one)

text: Assume that two r.v. $W$ and $Y|W=w$ with (1) $W \sim \text{N}(\mu_w,\sigma_w^2)$ (iid) (2) $Y|W=w \sim \text{N}(w,\sigma_y^2)$ (iid) Further we only observe $Y$ if $Y$ is less then $W$, i.e., (3) $Y|Y\le W$ **Goal:** Find the pdf of the censored observations, i.e., of $Y|Y\le W$ and from that deduce the uncensored pdf and the first two moments (so i.m.h.o. we have to find$f_Y(y)$). The first two moments of this uncensored pdf are supposed to depend upon $E(Y|Y\le W)$ and $Var(Y|Y\le W)$. * * * By definition of _conditional pdf_ we have that: (4) $f_{Y|W}(y|W = w)= \frac{f_{Y,W}(y,w)}{f_W(w)}$ Next, the definition of a _truncated density_ gives for a abitrary value of $W$: (5) $ f_{Y|Y\le W}(y|y\le w) = \frac{f_Y(y)}{P(Y\le W)}$ * * * I would simply rewrite (4) to $f_{Y|W}(y|W = w)f_W(w) = f_{Y,W}(y,w)$ then integration over $f_{Y,W}(y,w)$ w.r.t $w$ should yield $f_Y(y)$, i.e., (a) $\int_{-\infty}^{\infty} f_{Y,W}(y,w) dw = \int_{-\infty}^{\infty} f_Y(y|W = w)f_W(w) dw = f_Y(y)$ Plugin in $f_Y(y)$ into (5), ($P(Y\le W)$ will also be given by $f_Y(y)$) I will se how the moments of $f_{Y|Y\le W}(y|y\le w)$ will look and how the moments of $f_Y(y)$ depend upon them. So (a) will look like $f_Y(y) = \int_{-\infty}^{\infty}\frac{1}{\sqrt{2\pi\sigma^2_y}}\text{exp}\big(-\frac{(y-w)^2}{2\sigma_y^2}\big)\frac{1}{\sqrt{2\pi\sigma^2_w}}\text{exp}\big(-\frac{(w-\mu_w)^2}{2\sigma_w^2}\big)dw$ Except for the $w$ in the first $\text{exp}$, this looks very easy but since there is a $w$ im a little bit stuck how to solve this...


title: Verification of poisson approximation to hypergeometric distribution

text: How can I verify that $\lim_{N,M,K \to \infty, \frac{M}{N} \to 0, \frac{KM}{N} \to \lambda} \frac{\binom{M}{x}\binom{N-M}{K-x}}{\binom{N}{K}} = \frac{\lambda^x}{x!}e^{-\lambda}$, **without** using **Stirling's formula** or the **Poisson approximation to the Binomial**? I have been stuck on this problem for a while, because I don't know how to divide up the terms and factorials without using the help of prior results! Any help would be appreciated. Thanks in advance.


title: comparing simple(OLS) and multilevel regressions

text: Please, help me to understand how to compare two models - simple OLS and multilevel regression? Below you can see their pictures. In the simple regression there are 7 predictors including binary variable centr which is referred to central district of the city. In multilevel model I exclude centr and create the second level according to the district of the city (9 districts) Could I compare var(residual) in multilevel with residual MS in simple regression? Or maybe I can use AIC and BIC? Do I correctly understand that in simple regression we underestmate standard errors? Sorry for my chaotic question. It is my homework task and I am a bit confused with it. ![enter image description here](http://i.stack.imgur.com/p9oUZ.png) ![enter image description here](http://i.stack.imgur.com/ZaOo6.png)


title: Regression: Test a hypothesis about the product of beta and mean x

text: Two questions: a. Consider the regression y[i] = a + b x[i] + e[i]. Let bhat=(ols) estimate of b and xbar=sample mean of x. How to test the hypothesis (b*xu)=0? xu is the population mean of x. b. Consider the regression y[i] = a + b x[i] + c z[i] + e[i]. Let bhat and chat be the (ols) estimates of b and c, and xbar and zbar the sample means of x and z. How to test the hypothesis (b*xu + c*zu)=0? zu is the population mean of z. My first idea would simply be to set up a t-test: (bhat*xbar - 0) / SE(bhat*xbar). With SE(bhat*xbar) = V(bhat*xbar)^0.5. Is this appropriate? Are there better ways? Been looking for similar problems in a few statistical texts and on the web, but without luck. So in a. one would seem to need to work out V(bhat*xbar) and in b. V(bhat*xbar+chat*zbar). I think I can work out a) treating x as non-random. No responses, I note. Maybe the reason is that the question looks as if it is "self-study", but not tagged as such. (Self-study is here defined as "A routine question from a textbook, course, or test used for a class or self- study". In the info about the self-study tag, one can read "Be honest about the source of the question". The question does not fall under the definition as far as I can understand, therefore tagging it as self-study would be incorrect, maybe even dishonest. For example, I am not a student or enrolled in a course. If it is a routine question I would have expected to find something resembling it somewhere; but it might just be an inability to spot it. I am not even sure about what to search for exactly.) Anyway, I would most grateful if someone could direct me to a source where similar problems are discussed, in particular ones similar to b. If there is something fundamentally stupid about the questions, please point that out.


title: Matching results completely different in R and Stata

text: I am trying to match treatment with the control using the `R Matching package` for the data set nuclearplants (which is available in `optmatch package in R`). Following is the code for matching:               library (optmatch)# for the dataset nuclearplants     library(Matching) # for performing matching     data(nuclearplants)     # I am using Mahalanobis matching with covariates t1 and t2 and treatment pr     cov<-cbind(nuclearplants$t1,nuclearplants$t2)     sd<-Match(Y=NULL,Tr=nuclearplants$pr,X=cov, Weight.matrix = 2 )          matched<-rbind(nuclearplants[sd$index.treated,],nuclearplants[sd$index.control,])     head(matched)                      cost  date t1 t2  cap pr ne ct bw cum.n pt     A 443.22 67.33 10 85 1065  1  0  1  0     1  0     B 642.23 68.00 11 78 1065  1  1  1  0    12  0     C 457.12 68.42 15 55  822  1  0  0  0     5  0     D 289.66 68.42 15 76  530  1  0  1  0     2  0     E 490.88 68.92 16 59 1050  1  0  0  0     8  0     F 665.99 70.92 22 57  828  1  1  0  0    20  0     > sd$index.treated           [1]  3  5  9 18 20 22 24 30 31 32         > sd$index.control      [1]  2 10  8 19 23 26 26 10 10  2      Now, I am trying to see whether the `stata` also gives the same results for the same data and the same methodology. For this purpose, I use `psmatch 2`               psmatch2 pr, mahalanobis(t1 t2)      I was very shocked to see that for the same approach, these two programs give different results. Below is the summary results for the stata:               tre stata control     3   7     5   7     9   6     18  15     20  17     22  19     24  14     30  7     31  7     32  2          tr: treatment observations (row number), stata control=control observations for Stata (row number)      I was wondering whether I missed something implicit in the code.


title: SPSS: Determining the direction

text: Background: Group A the control group uses a website 1. Group B the treatment group uses website B. The expected result is that Group B will have a higher understanding than group A. Hypothesis: The treatment group will show increased levels of understanding. A 1(Strongly Agree) - 7(Strongly disagree) likert survey was used, to measure understanding after the groups used either website 1 or 2. I know using a two sample t-test will allow me to determine whether a significant difference was present. But how do I determine the direction of the difference through using SPSS? Still new to this. Thanks for the help.


title: How to assess relationship between responses to multiple likert items between two raters when there are multiple rater pairs?

text: I have two data sets in which the same Likert scale questions were answered by two different responders. I want to calculate which level they are correlated. My data set is like               WorkerID    Question    Worker_answer     Employer_answer_about_his_worker      1               1              3                       2     1               2              4                       3     1               3              1                       5     1               4              2                       4     .               .              .                       .     1               15             3                       2     2               1              4                       2     2               2              5                       5     .               .              .                       .     2               15             2                       3


I have calulated an OLS with and without a constant term. However, besides that the values are different I haven`t really found anything valueable? Therefore my question is: Whats the biggest difference in calculating a simple regression model with/without a constant term? I really appreciate your reply!


title: What are the main differences between K-means and K-nearest neighbours?

text: I know that k-means is unsupervised and is used for clustering etc and that k-NN is supervised. But I wanted to know concrete differences between the two?


title: Monotonic and non-monotonic patterns of missing values: how do they look like?

text: I was reading the user's guide of SPSS on missing value analysis and I found the terms monotonic and non-monotonic pattern of missing values. The terms are not quite clear to me. They used `telco_missing.sav` data file to show a graph from which they comment: > If the data are monotone, then all missing cells and nonmissing cells in the > chart will be contiguous; that is, there will be no “islands” of nonmissing > cells in the lower right portion of the chart and no “islands” of missing > cells in the upper left portion of the chart. I didn't either understand that. I hope someone will be kind enough to show me how a monotone and non-monotone patterns does look like and tell me how do I understand it a little more clearly.


title: Error in JAGS code

text: I'm trying to run the following code in JAGS               model {     for (i in 1:n){         z[i] ~ dnorm(X[i], tau)         X[i] <- D1*x1[i] + D2*x2[i] + D3*x3[i] + D4*x4[i]     }          denom <- 1 + sum(exp(p))     D1 <- (exp(p[1]))/denom     D2 <- (exp(p[2]))/denom     D3 <- (exp(p[3]))/denom     D4 <- 1/denom          # Prior for Fixed Effects     p[1] ~ dnorm(0,.01)     p[2] ~ dnorm(0,.01)     p[3] ~ dnorm(0,.01)          # Prior for precision of Delta     tau2 <- pow(sigma2, -2)     sigma2 ~ dunif(0, 100)          # Prior for Variance of Z     tau <- pow(sigma, -2)     sigma ~ dunif(0, 100)          }      In R I run               # Data     z <- d1$transfers.cap         x1 <- ifelse(d1$LM==1, 1, 0)     x2 <- ifelse(d1$LM==2, 1, 0)         x3 <- ifelse(d1$LM==3, 1, 0)     x4 <- ifelse(d1$LM==4, 1, 0)     n <- dim(d1)[1]          ex.data <- list(z=z, x1=x1, x2=x2, x3=x3, x4=x4, n=n)     out <- jags.model(file="SEM.bug", data=ex.data, n.chains=5, n.adapt=1000)      But I get the error message               Compiling model graph        Resolving undeclared variables        Allocating nodes     Deleting model          Error in jags.model(file = "SEM.bug", data = ex.data, n.chains = 5, n.adapt = 1000) : RUNTIME ERROR:     Invalid vector argument to exp      How can I solve this?


title: Homogenity of Variances - F-test and Levene Test yield different results - which one to trust?

text: I want to perform a t-test on two independent samples with each n = 50. to check the assumptions, I ran an F-test for homogenity, which is significant. then I used SPSS which does a levene test by default and this levene test was not significant. I don't have the data here, but n = 50 (data are rental prices) is probably normal distributed. which test should I trust and why? thanks in advance for your answers!


title: Resources or methods for determining if a technique is in current use?

text: I know that this question is hard to answer, but I was wondering how you would determine if a statistical technique was in current use, or if it had a heyday a couple of decades ago but it's either been superseded by newer techniques or perhaps it was even proven to be inadequate and abandoned? For example, I was looking at a time series with R and stumbled onto the EMD package (Empirical Mode Decomposition). The technique looks like interesting, but as I google around, it seems that it never really became popular. Was it confined to the neural research community? Did wavelets prove to be superior? Did it prove to be a bit hard to interpret? Perhaps there can be no definitive answer, but perhaps there are signs you look for?


title: Assessing difference observations of two samples

text: I have two continous samples: background and data. I want to check whether data observations are significantly greater than background. I used U-Test to calculate p-value, but now I want to find a pair of background/data with the biggest `background >> data` between them. In practice I have a formula to which calculates some stuff and return one value per observation. I want to show that if I apply this formula the result is significantly bigger for relevant data (non-background) and then among 6 formulas I want to choose the one which gives a biggest difference in results.


title: How to encode categorical variables for neural networks

text: In regression, you encode a categorical variable with n possible values using n-1 indicator variables. How about for neural networks?


title: Create 3 correlated variables with pre-specified correlations

text: I start with a vector of random numbers sampled from a normal distribution:               R<-rnorm(100, mean=0, sd=30)      I would now like to create 3 variables that are correlated with each other with a pre-specified correlation. In addition I would like to have these three variables correlated with `R` with a pre-specified correlation. For example `A`, `B`, `C` would have correlation 0.7 with each other, and `A`, `B`, `C` would have correlation 0.6 with `R`. I.e. I am looking for the following covariance matrix:                   R   A   B   C     R   1  .6  .6  .6     A  .6   1  .7  .7     B  .6  .7   1  .7     C  .6  .7  .7   1      How can this be done in R?


title: How important is it to transform variable for Cox Proportional Hazards?

text: I am building a Cox PH model for recurrence-free survival from breast cancer. One of my variables, the number of involved lymph nodes, is highly skewed with many more "0" values than positive integers. It seems to me that this would make any transformation ineffective because this variable cannot be negative. Here is a stem plot of the lymph node data:                0 | 00000000000000000000000000000000000000000000000000000000000000000000+195      2 | 00000000000000000000000000000000000000000000000000      4 | 00000000000000000000000000000      6 | 000000000000000      8 | 000000000000     10 | 0000000000     12 | 00000     14 | 0000     16 | 0     18 |      20 | 0     22 |      24 | 00     26 | 0     28 | 0     30 |      32 | 0      A log(X+1) transformation doesn't do anything to normalize the data because there is no data less than 0. How important is it to transform this variable? Would using only values <20 and log-transforming them help more? This must be a common problem; are there any creative solutions to it? Thanks for your help!


title: Fast diagonalization for coordinate descent

text: IN order to perform efficient coordinate descent, I am using SVD to diagonalize my system $$y=Ax$$ to $$y=USV^Tx\\ U^Ty=y_{new}=S(V^Tx)=Sx_{new}\\ y_{new}=Sx_{new}$$ where $S$ is the diagonal matrix of singular values. Now I can solve the system in parallel by just $x_i=\frac{y_i}{s_i}$. The problem is, that SVD has the complexity of $O(N^3)$ and I am looking for a more efficient process. However, maybe I am completely wrong here and a less complex preconditioning process can be used for coordinate descent? Thanks a lot for any hint!


title: Correlation between two variables measured on a “strongly agree” to “strongly disagree” scale

text: The questions on a survey asked:   * Do you actively participate in a study group?   * Do you think the class is going too quickly? For both, the responses are one of the following: strongly agree, agree, neutral, disagree, strongly disagree So I want to analyze if there is a correlation between those who are in a study group and those who think the class is going to quickly. So I have two columns in a data frame. The first column is labeled "group" and the second is "fast". I have converted the group variable into a numerical variable, like so: strongly agree = 5, agree = 4, neutral = 3, disagree = 2, strongly disagree = 1 So now I have a data frame with two columns, one full of numbers and the other still full of the original answers ("strongly agree", "agree", etc). I have found the means of every option quickly, now I just need to see if there is a statistical significance, but I am clueless. How exactly should I calculate the p-value on this? I have tried several methods but the p-value seems wrong. Sorry if this is easy stuff, I think I have made this way more complicated in my mind than it should be and I appreciate any help.               quickly <- CSExperiencesAllWithHeaders$CEQuickly         groups <-CSExperiencesAllWithHeaders$CEStudyGroup          levels(groups) <- (c(levels(groups), 5, 4, 3, 2, 1))     groups[groups == "strongly agree"] <- 5     groups[groups == "agree"] <- 4     groups[groups == "neutral"] <- 3     groups[groups == "disagree"] <-2     groups[groups == "strongly disagree"] <- 1     groups[groups == ""] <- NA     groups[groups == "N/A"] <- NA     quickly[quickly == "N/A"] <- NA     quickly[quickly == ""] <- NA          groups <- factor(groups)     quickly <- factor(quickly)     analysis3 <- data.frame(groups,quickly)     analysis3 <- na.omit(analysis3)     analysis3$groups <- as.numeric(as.character(analysis3$groups))          sagree2 <- subset(analysis3, quickly  == "strongly agree")     agree2 <- subset(analysis3, quickly == "agree")     neutral2 <- subset(analysis3, quickly == "neutral")     disagree2 <- subset(analysis3, quickly == "disagree")     sdisagree2 <- subset(analysis3, quickly == "strongly disagree")          meansagree2 <- mean(sagree2$groups)         meanagree2 <- mean(agree2$groups)     meanneutral2 <- mean(neutral2$groups)         meandisagree2 <- mean(disagree2$groups)     meansdisagree2 <- mean(sdisagree2$groups)          barplot(c(meansagree2, meanagree2, meanneutral2, meandisagree2,                meansdisagree2),             main = "Those Who Think Class is Too Quick: In Study Groups?",             names.arg=c("Strongly Agree","Agree","Neutral","Disagree",                          "Strongly Disagree"),             xlab = "Class too Quick?",             ylab = "In a Study Group?")      all this code creates this data frame (I only took the top of the data frame since the real one is over 1000 columns):                   groups  quickly     1   5   'strongly disagree'     2   4   'strongly agree'     3   1   'disagree'     4   1   'disagree'     5   4   'strongly disagree'     6   2   'strongly disagree'     7   1   'neutral'     8   2   'disagree'     9   1   'strongly disagree'     10  2   'strongly disagree'     11  1   'strongly disagree'     12  2   'neutral'     13  5   'disagree'     14  2   'disagree'     15  4   'neutral'     16  2   'disagree'     17  5   'disagree'     18  5   'neutral'     19  4   'strongly disagree'     20  2   'strongly disagree'     21  3   'disagree'     22  1   'strongly disagree'     23  4   'strongly agree'     24  1   'strongly disagree'     26  5   'strongly disagree'     27  1   'strongly disagree'     28  5   'disagree'     29  5   'agree'      This is what I get when I use the dput function: structure(list(groups = c(5, 4, 1, 1, 4, 2, 1, 2, 1, 2, 1, 2, 5, 2, 4, 2, 5, 5, 4, 2, 3, 1, 4, 1, 5, 1, 5, 5, 5, 5), quickly = structure(c(5L, 4L, 2L, 2L, 5L, 5L, 3L, 2L, 5L, 5L, 5L, 3L, 2L, 2L, 3L, 2L, 2L, 3L, 5L, 5L, 2L, 5L, 4L, 5L, 5L, 5L, 2L, 1L, 2L, 3L), .Label = c("agree", "disagree", "neutral", "strongly agree", "strongly disagree"), class = "factor"), qui_fact = structure(c(5L, 1L, 4L, 4L, 5L, 5L, 3L, 4L, 5L, 5L, 5L, 3L, 4L, 4L, 3L, 4L, 4L, 3L, 5L, 5L, 4L, 5L, 1L, 5L, 5L, 5L, 4L, 2L, 4L, 3L), .Label = c("strongly agree", "agree", "neutral", "disagree", "strongly disagree"), class = "factor"), qui_num = c(5, 1, 4, 4, 5, 5, 3, 4, 5, 5, 5, 3, 4, 4, 3, 4, 4, 3, 5, 5, 4, 5, 1, 5, 5, 5, 4, 2, 4, 3)), .Names = c("groups", "quickly", "qui_fact", "qui_num"), na.action = structure(c(25L, 31L, 37L, 38L, 86L, 91L, 148L, 209L, 270L, 280L, 285L, 328L, 338L, 340L, 410L, 424L, 456L, 460L, 461L, 480L, 568L, 587L, 593L, 596L, 599L, 600L, 607L, 621L, 658L, 700L, 717L, 731L, 758L, 776L, 827L, 837L, 849L, 862L, 864L, 896L, 899L, 909L, 921L, 946L, 963L, 966L, 977L, 994L, 1007L, 1012L, 1074L, 1079L), .Names = c("25", "31", "37", "38", "86", "91", "148", "209", "270", "280", "285", "328", "338", "340", "410", "424", "456", "460", "461", "480", "568", "587", "593", "596", "599", "600", "607", "621", "658", "700", "717", "731", "758", "776", "827", "837", "849", "862", "864", "896", "899", "909", "921", "946", "963", "966", "977", "994", "1007", "1012", "1074", "1079"), class = "omit"), row.names = c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L, 23L, 24L, 26L, 27L, 28L, 29L, 30L, 32L), class = "data.frame")


title: Theorem about ARMA process existing

text: I have met the theorem: "For each stationary time series there is ARMA model with finite number of parameters." Could you indicate official source with proof (bibligraphy) of this theorem?


title: Controlling for multiple, nominal variables in a linear model

text: Here's my experiment: I'm testing the physiological effects of a treatment on a mice (measuring continuous variables such as total mass, fat mass, bone density, etc). I want to test the significance of this treatment on weight gain parameters. What I've done so far: in R, I've fitted a logistic glm model, using the binary factor indicating treatment/control group as the response (I eventually want to see what weight measurements are predictive of treatment group- membership). However, I also have repeated measures at (for now) 2 time points (~3 weeks apart) and a sex variable as well. These are (uninteresting) factors effecting weight and weight gain. My question is, how do I control for factors that I know are significantly effecting weight gain? The problem seems to be that I want to use continuous measurements to predict a classifier (treatment group) but this is probably getting swamped out by trivial attributes that are also effecting these continuous measurements (age and sex). Other than just splitting up the data into subgroups, is there are way to deal with this sort of 'mixed effect' model? A reference would be greatly appreciated. Thanks!


title: Bayesian inference of parameters: residuals are independent but not normally distributed

text: I would like to compute belief intervals (confidence intervals; CI) for the parameters of an environmental dynamic model within the Bayes' theorem. The measurement model of the data is $$ d(t)=y(\theta,t)+\epsilon, $$ where $\epsilon$~Norm(0, $\sigma$) is the assumed measurement error, $y$ is the dynamic system output dependent on parameter $\theta$ and $d(t)$ is the observed data at time t. Because the measurment error is assumed normally distributed, the likelihood function is a Gaussian distribution. The prior for $\theta$ is non-informative (i.e. uniform distribution). The aim is to compare CI estimated from Bayes’ formula and the classical linear approximation (covariance matrix estimated from the Hessian inversion): I would like to maintain the “least-square” objective function in order to perform this comparison. The measurement model assumes that errors are i.i.d. (independent and identically distributed). What is the meaning of the CI if the residuals are not normally distributed but are still "independent" (without auto- correlation)? In my case the distribution of residuals has long tails (e.g. Laplace distribution). In the reference book of Jaynes et al. (2003), pp.592–593, it is said that: _Measurement errors in physical experiments are often modeled by a normal distribution. This use of a normal distribution does not imply that one is assuming the measurement errors are normally distributed; rather using the normal distribution produces the most conservative **predictions** possible given only knowledge about the mean and variance of the errors._ I understand that Jaynes' statement refers to the prediction uncertainty. What about CI of $\theta$, apart of $\sigma$? Should $CI_{\theta}$ be considered overestimated and conservative? Or the opposite since we do not consider a statistical model for residuals with long tails (i.e. the likelihood is proportional to the posterior)? Is the CI estimation totally unreliable because normality of residuals is not guaranteed or is it "unreliable" only to the “degree” the normality assumption is violated?


title: Business forecasting

text: Data can have trends, cycles, seasonal pattern repeats and irregular components. How do these components impact on the ordinary multiple regression with dummy variables?


title: Implementation of stratified crossvalidation for classification problem

text: Is this all I have to do to implement 2 fold stratified cross validation:   1. Assign random number from uniform distribution [0...1] to each row.   2. Assign row to fold 1 if random number <= 0.5 otherwise assign to fold 2 Does this ensure that the class priors are roughly the same in each fold?


title: pooled regression in eviews

text: Hey I have to regress abnormal returns for 202 companies over a period of 10 years (between 2004-2013) against their respective PE ratios and a Profitability measure for the same time period. Thus I have to run a pooled regression. I have read online on how to do it...but still am not sure of the process and also did not find any appropriate examples, or one similar to my case. Here my dependent variable are the returns of the 202 stocks over a 10 yr time period. And the two independent/ explanatory variables are their PE ratios, current and next year and profitability measures respectively (both of which are also considered fr the same time period of 10 years). So please could anyone tell me whether it is possible to run this regression in eviews? I was confused regrading what should i select in the estimation box options. Firstly, what should I consider as my common coefficient/cross section coefficient..in the estimation box options. Also should I go for fixed effects/common/ random effects intercept option? The equation I constructed to test whether there are any PE, Expected PE or time lag PE effects on Abnormal Returns is as follows: Abnormal Return (t) = α+β1 (PE t-1)+β2(PE t)+β3 (PE t+1)+ β4(Profitability t)+ε Please help me out anyone by telling me how can I run a pooled regression for my dataset and equation..or is it even possible..to take these data from three different sheets and accordingly run it because otherwise I have to take mean returns of all companies and mean PE of all companies and run the regression which would be less satisfactory to my dissertation research. Thanks alot in advance. Regards, Samiha


title: Expected value of a series of random variables in a markov chain

text: I have a Markov Chain such that $X_n = max(X_{n-1}+\xi _n,0)$ where the $\xi_n$ series is independent and identically distributed. I want to show that if $\mathbb E(\xi_n) > 0$ (where $\mathbb E(\xi_n)$ is the expected value of $\xi_n$) then$\frac{X_n}{n}$ tends to $\mathbb E(\xi_n)$ as n approaches infinity for any choice of $X_0$. And... I have no idea how to show this. I know that repeated application of the above will yield $X_n = max(X_0+\sum\limits_{i=1}^n\xi_i, \sum\limits_{i=2}^n\xi_i, \sum\limits_{i=3}^n\xi_i, ..., 0)$, but I'm stuck here.


title: Computing Confidence in A/B testing

text: First of all, I'm fairly new to Decision Theory and my question might be rather too simple for the experts. I am having some issues in computing confidence for A/B Testing. I will first explain the scenario and then the technique to compute confidence. It's a typical scenario extending from clinical trials, where I test two headlines for an article to see which one's doing better. Let's say that I have two headlines and I do the A/B test for 20 minutes and see which one is better and by how much confidence. I compute the CTRs of the headlines (= clicks/views) after 20 minutes and see if there's a "significant" difference between the two. I use the Anscombe quartet for that and if the difference in CTRs is significant enough, I pick a winning headline. To compute the confidence, I model the CTRs of the headlines as beta distributions with $\alpha = ($VIEWS$_a \times$ CTR$_a)$ and $\beta = ($VIEWS$_a \times (1 - $CTR$_a))$, which obviously assumes that the prior is a uniform distribution. Once I have the two beta distributions, I compute the confidence using Monte Carlo Sampling of the two distributions. I set my cut- off as 80% for the overlap to be able to "confidently" tell that one headline is better than the other. However, since the kurtosis of the beta distribution increases non-linearly, for a given CTR, it increases non-linearly with the number of views. So, for very small differences in CTRs, I can end up with very high confidence if the number of views is very high for a given article. Is there anything I can do to fix this? Use some other distribution with some other confidence measure? I also want the distribution to make good intuitive sense and I'm not really sure if Beta is the best that I can choose for this. Thanks in advance!


title: Bootstrap on correlation coefficients in R

text: For each round of bootstrapping (say $n=1000$), I would like to select the genotypes (the columns) at random without replacement from the original dataset. With this I would like to form two groups with 10 and 10 genotypes each to match the genotypes in $X$ and $Y$ subsets respectively. With the original dataset (correlation between $X$ and $Y$) and the bootstrapped datasets (correlation between $X'$ and $Y'$) draw a histogram as follows to see how many percent of the 1000 pairs derived from random permuting exhibit lower correlations than the original dataset. A sample histogram is shown below. ![enter image description here](http://i.stack.imgur.com/YE06t.png) Maybe for such a sample data                set.seed(1)        x <- matrix(rnorm(1000), nrow=100, ncol=10)       y <- matrix(rnorm(1000), nrow=100, ncol=10).      I found a useful MATLAB code which I have added here to bootstrap correlation coefficients. Is there something similar in R to obtain a similar plot?


title: Power of Seasonal Kendall Trend test

text: I was asked about doing a power analysis of a seasonal kendall trend test. I feel like that would be really difficult to do and I haven't found any documentation or software on how to do it. Is there some way to estimate the required sample size using regression? Hopefully someone can point me in the right direction. Thanks!


title: "Normalized" standard deviation

text: I posted a question earlier but failed miserably in trying to explain what I am looking for (thanks to those who tried to help me anyway). Will try again, starting with a clean sheet. Standard deviations are sensitive to scale. Since I am trying to perform a statistical test where the best result is predicted by the lowest standard deviation amongst different data sets, is there a way to "normalize" it for scale, or use a different standard-deviation-type test altogether? Unfortunately dividing the resulting standard deviation by the mean in my case does not work, as the mean is almost always close to zero. Thanks


title: Prediction interval on X with Poisson

text: What I'm trying to do is predict some biodata based on time. Because the biodata are not normally distributed, I have modeled them with a poisson using a log link function ( E[y] = exp[BX]). What I'm trying to do is estimate the time at which the predicted curve reaches a certain threshold. (i.e., y=1 when x=?). Obtaining a point estimate is fairly straightforward, but I would like to place a confidence interval around that time (e.g., x = 5+-1). My initial idea was to flip the x and the y and generate prediction intervals as one normally would. However, x is normally distributed while y is poisson, so that's not a good solution. Any other ideas?


title: Confidence of uniformity from a sample of a binary state collection

text: The system:   1. I have a bag filled with 1,000 balls   2. All of them are either red or blue I want to gain some amount of confidence that all the balls are blue, but I'm very lazy and I only want to examine a sample as small as possible. **Question: How large must the sample be to say with P confidence that N balls are blue?** I'm looking for a generic answer, but here are two specific questions:   1. How large must the sample be to say with 95% confidence that all 1,000 of the balls are blue?   2. How large must the sample be to say with 95% confidence that 990 (99%) of the 1,000 balls are blue? It's been a while with no response, so if my question is missing information or just bad for some reason, please let me know how I can improve it.


title: Threshold for Fisher linear classifier

text: The fisher linear classifier for two classes is a classifier with this discriminant function: $h(x) = V^{T}X + v_0$ where $V = \left[ \frac{1}{2}\Sigma_1 + \frac{1}{2}\Sigma_2\right]^{-1}(M_2-M_1)$ and $M_1$, $M_2$ are means and $\Sigma_1$,$\Sigma_2$ are covariances of the classes. $V$ can be calculated easily but the fisher criterion cannot give us the optimum $v_0$. Then, how we can find the optimum $v_0$ analytically if we are using fisher classifier?


title: Non-parametric test for unimodality of 2-d dataset

text: I'm looking for a method to determine if a 2-d dataset forms a single cluster or is split into two/multiple clusters - and what's the likelihood or p-value for this. Isolated data should not be considered as different cluster. Also, data has a non-parametric distribution. My approach is to apply an unimodal / bimodal test. However, I didn't find any modality statistical test for 2-d non-parametric data. An alternative approach is to check the modality of the projections of the data on the two axis (this becomes a 1-d modality test). If the distributions of the projections are unimodal on both axis, then the 2-d dataset has a unimodal distribution. If this alternative approach is correct, then my question becomes: what unimodality test to use for 1-d non-parametric data? - I already tried the DipTest R implementation, but it doesn't work correctly in the non-parametric case (see the example below). An example illustrating the alternative approach below: (first two plots - histograms of the projection of the data on the OX and OY axes; the third plot - scatterplot of the 2-d dataset)   * Example for unimodal distribution - unimodal distribution on both axes (inferred from the histograms) !! inconsistent results from DipTest: pval_x = 0.06 & pval_y = 0.05 (which means significant bimodlity on both axis - which is not true) ![unimodal](http://swarm.cs.pub.ro/~mdragan/gsoc2013/genevalidator/stat/unimodal.png)   * Example for multi-modal distribution  ![bimodal](http://swarm.cs.pub.ro/~mdragan/gsoc2013/genevalidator/stat/bimodal.png) Thanks very much for your help, Monica


title: Calculating the error of Bayes classifier analytically

text: If two classes $w_1$ and $w_2$ have normal distribution with known parameters ($M_1$, $M_2$ as their means and $\Sigma_1$,$\Sigma_2$ are their covariances) how we can calculate error of the Bayes classifier for them theorically? Also suppose the variables are in N-dimensional space. **Note:** A copy of this question is also available at http://math.stackexchange.com/q/11891/4051 that is still unanswered. If any of these question get answered, the other one will be deleted.


title: Generated regressor and interaction term

text: I want to estimate the following model $y_{it} = \gamma_{0} + \gamma_{1} x_{it} + \gamma_{2} \theta_{i} x_{it} + u_{it}$ where the $\theta_i$ is unknown but I can estimate it in another regression. Plugging in that estimate and rewriting the equation, I get $y_{it} = \gamma_{0} + \gamma_{1} x_{it} + \gamma_{2} \hat \theta_{i} x_{it} + u_{it} + \gamma_{2} x_{it} (\theta_i - \hat \theta_{i})$ The last two terms are obviously unobserved, and the last one correlates with $x_{it}$, so my estimates are inconsistent. I am hoping that this is a common problem that has some established fix. Does anyone know literature on that? Thank you in advance.


title: Is there a package or software in R for calculating regression by the York method?

text: I have elemental concentrations for 18 elements in 36 samples, determined by neutron activation analysis. I would like to calculate correlation coefficients, etc. by regression analysis according to the York method.


title: Applying post hoc paired t-tests with FDR correction when you have two significant main effects but no interaction - is this permitted?

text: I'd be grateful for your help. I have run a repeated measures ANOVA in SPSS. I have two factors - I and II. Factor I has 2 levels (a) and (b) and factor II has 3 levels (c) (d) and (e). I have identified significant Main effects of Factor I and Factor II. However, there is no significant interaction between Factor I and Factor II. Am I allowed to do the following post hoc paired t-tests (*=significant) in the absence of an interaction but when there are two main effects as per the ANOVA? Factor I, Level a/Factor II, Level c v Factor I, Level a/Factor II, Level d * Factor I, Level a/Factor II, Level c v Factor I, Level a/Factor II, Level e * Factor I, Level b/Factor II, Level c v Factor I, Level b/Factor II, Level d ns Factor I, Level b/Factor II, Level c v Factor I, Level b/Factor II, Level e * Thanks for your help!


title: Wald-Wolfowitz Runs Test for Normality Assumptions Testing

text: Some guys told me that it's appropriate to use Wald-Wolfowitz Runs Test as a normality test (like Shapiro-Wilk's or Kolmogorov-Smirnov...). Do you think this is good way to test normality assumptions?


title: Confidence Interval for Percentile for Empirical Distribution

text: I have a bunch of raw data values that are dollar amounts and I want to find a confidence interval for a percentile of that data. Is there a formula for such a confidence interval? Thanks for any help


title: When is a second hidden layer needed in feed forward neural networks?

text: I'm using a feed forward neural network to approximate a function with 24 inputs, and 3 outputs. Most of the literature suggests that a single layer neural network with a sufficient number of hidden neurons will provide a good approximation for most problems, and that adding a second or third layer yields little benefit. However I have optimized a single layer, and a multi-layer neural network and my multi-layer network is much better. For the single layer network I performed a sweep of 1 to 80 neurons, retraining the network each time, and plotting the performance. After about 30 neurons the performance converged. For the multi-layer network I used the genetic algorithm to select the number of neurons in the first and second layer. This resulted in a much better performance. What I would like to know is why this happened, considering that most of the literature suggests that 1 layer is enough. Is there a particular type of problem that requires more than one layer? Does this suggest that the function being approximated is discontinuous, not well-defined, jagged (not smooth), or all/ a mix of the above? Or does it suggest something else, or nothing at all? I know that when used for classification a multi-layer neural network can classify data that is not linearly separable, but I'm using the network for function approximation.


title: Log likelihoods of subsamples don't sum up to full sample log-likelihood

text: When splitting the sample into segments, I would expect the likelihood of the segments to sum up to at least the likelihood of the full sample (constraining the parameters to be the same). Or did I misunderstand? The code example is in `R`, using the `rugarch` package. The `rmgarch` package is only used to get the data. First I calculated the log-likelihood of a GARCH model that was fitted to the full sample. Then I split the sample in two parts and and calculated the log- likelihood for each segment while fixing the parameters to be the same as the ones obtained by step 1. I would expect the subsample log likelihoods to sum up to the log likelihood of the full sample. But this doesn't happen. When the segments-GARCH models are not filtered with fixed parameters, but instead allowed to estimate their own parameter estimates, their sum is still below the likelihood of the full sample. Actually it should be equal or bigger than the full likelihood. Did I misunderstand the concept of ML??? Best, Jo               require('rugarch')     require('rmgarch')     data(dji30retw)     dat <- dji30retw$AA          ##### 1: get log-likelihood and parameters of full-sample model     spec1 <- ugarchspec(          variance.model = list(model = "sGARCH", garchOrder = c(1,1), submodel = NULL)          , mean.model = list( armaOrder =c(0,0) ) , distribution.model = "norm" )     garch1 <- ugarchfit(spec=spec1, data=dat )     l1 = likelihood(garch1)          spec2 = spec1     setfixed(spec2)<-as.list(coef(garch1))          ##### 2: filter two sample-segments individually using the parameters from step 1     garch2 = ugarchfilter(spec = spec2, data = dat[1:344])     l2 = likelihood(garch2)          garch3 <- ugarchfilter(spec=spec2, data=dat[345:length(dat)]  )     l3 = likelihood(garch3)          l1        # 1901.706     l2+l3     # 1898.931     # should be the same in theory          ##### 3: now fit the individual segments     garch5 <- ugarchfit(spec=spec1, data=dat[1:344] )     l5 = likelihood(garch5)          garch6 <- ugarchfit(spec=spec1, data=dat[345:length(dat)] )     l6 = likelihood(garch6)          l1        # 1901.706     l5+l6     # 1899.254     # the sum of the individual log-likelihoods shouldn't be smaller than the full-sample log-likelihood!


title: Computing best subset of predictors for linear regression

text: For the selection of predictors in multivariate linear regression with $p$ suitable predictors, what methods are available to find an 'optimal' subset of the predictors without explicitly testing all $2^p$ subsets? In 'Applied Survival Analysis,' Hosmer & Lemeshow make reference to Kuk's method, but I cannot find the original paper. Can anyone describe this method, or, even better, a more modern technique? One may assume normally distributed errors.


title: Assessing clustering statistical significance

text: I read a good tutorial on the p-value. Everything there is clear for me, except how I can define the null hypothesis. I found an example (on page 6) in this paper using p-value to judge the statistical significance of the cutpoint he defined in his hierarchical clustering. The author proposed a permutation based approach: He permutes the data, evaluates the statistic of interest. He repeats this operation several times and obtains a distribution on the statistic from which he computes a p-value. I did a clustering using partitioning around medoids after I divided the problem using a canopy approach. I want to assess the statistical significance of the clusters I obtained by computing a p-value because, for some of these clusters, most of the features have 0 value on every sample. How can I define my null hypothesis for this purpose?


title: How to include $x$ and $x^2$ into regression, and whether to center them?

text: I want to include the term $x$ and its square $x^2$ (predictor variables) into a regression because I assume that low values of $x$ have a positive effect on the dependent variable and high values have a negative effect. The $x^2$ should capture the effect of the higher values. I therefore expect that the coefficient of $x$ will be positive and the coefficient of $x^2$ will be negative. Besides $x$, I also include other predictor variables. I read in some posts here that it is a good idea to center the variables in this case to avoid multicollinearity. When should you center your data & when should you standardize?   1. Should I center both variables seperately (at the mean) or should I only center $x$ and then take the square or should I only center $x^2$ and include the original $x$?   2. Is it a problem if $x$ is a count variable? In order to avoid $x$ being a count variable, I thought about dividing it by a theoretically defined area, for example 5 square kilometers. This should be a little bit similar to a point density calculation. However, I am afraid that in this situation my initial assumption about the sign of the coefficients would not hold anymore, as when $x=2$ and $x²=4$ $x= 2 / 5 \text{ km}^2$ = $0.4 \text{ km}^2$ but $x^2$ would then be smaller because $x^2= (2/5)^2= 0.16$.


title: Flowcharts to help selecting the proper analysis technique and test?

text: As someone who needs statistical knowledge but is not a formally trained statistician, I'd find it helpful to have a flowchart (or some kind of decision tree) to help me choose the correct approach to solve a particular problem (eg. "do you need this and know that and that and consider data to be normally distributed? Use technique X. If data is not normal, use Y or Z"). After some googling, I've seen several attempts of various coverage and quality (some not available at the moment). I've also seen similar flowcharts in statistics textbooks I've consulted in libraries. A bonus would be an interactive site that, besides just having a chart, would provide extra info (such as assumptions) and would point to how to perform those techniques in popular stat packages. "Need to do ANOVA in R? You need package X and here's a tutorial". I'm asking as a community wiki question in the hope there are better resources I couldn't find. Since statistics is a large subject, I think such a flowchart would be suitable for techniques that can be approached by someone who has beginner or intermediate-level knowledge. Anything more complicated would need someone with formal training.


title: Problem in evaluating naive Bayes

text: I am trying to understand naive Bayes and its application to text classification. I have a doubt or this may be my misconception. Suppose we have two categories "News" and "Sports" in which we need to classify any given document. Let the dictionary contain only 3 keywords ${news,football,tennis}$ with the following parameters \begin{align*} P(news/News)&=0.99,P(news/Sports)=0.01,P(tennis/Sports)=0.9,\\ P(tennis/News)&=0.1,P(football/Sports)=0.9,P(football/News)=0.1, \\ P(Sports)&=0.5 ,P(News)=0.5,P((news,football,tennis))=k; \end{align*} We get a document which has all the three keywords. So when we evaluate \begin{align} P(News/(news,football,tennis))=0.99\cdot0.1\cdot0.1\cdot0.5/k=0.00495/k\\ P(Sports/(news,football,tennis))=0.01\cdot0.9\cdot0.9\cdot0.5/k=0.00405/k \end{align} So the document is classified to "News" category, but intuitively we know that it should belong to "Sports" category.


title: A matching problem for 3 decks of cards

text: a recent posting on math.stackexchange reminds me of a somewhat less ambitious$^1$ question i have been meaning to ask. suppose 3 people each have a deck of $M$ cards numbered from 1 to $M$, in random order. together, they turn over one card at a time. on any such trial, a match occurs if at least 2 cards show the same number. let $X$ be the number of matches that occur in the $M$ trials. what is the probability distribution of $X$? [for starters, what is ${\rm P}(X=0$)?] an historical note: the solution for 2 players is well-known; _cf_ $ $ Feller, vol.1,3$^{rd}$ Ed., _p_ 107: $\quad {\rm P}(X=0) := p_{0} = 1 - 1 + \frac{1}{2!} - \frac{1}{3!} + \cdots \pm \frac{1}{M!}$, which, for $M$ not too small, is about 1/e. [for $M$ = 10, the two numbers already agree, on rounding, to 5 decimal places.] more generally, $\quad {\rm P}(X=k) := p_{k} = \frac{1}{k!} \{ 1 - 1 + \frac{1}{2!} - \frac{1}{3!} + \cdots \pm \frac{1}{(N-k)!}\} \sim \frac{e^{-1}}{k!}.$ [so as $M \to\infty, X$ tends in law to a poisson variable with mean 1.] for $k=M-1$ this gives: $p_{M-1} = 0$ [which is obvious]. also, $p_M= \frac{1}{M!}$. Feller obtains these formulae using the inclusion-exclusion principle. * * * $^1$ in the sense that it involves only 3 people, rather than everyone in the world, and they shuffle their decks only once, rather than repeatedly since the beginning of time.


title: Interpreting 2D correspondence analysis plots (Part II)

text: I'd like to ensure that I understand the process correctly. This is a follow- up question to http://stats.stackexchange.com/questions/3270/interpreting-2d-correspondence- analysis-plots               library(reshape)      library(ca)          df <- read.csv(file="http://www.bertelsen.ca/R/smokers.csv")     colnames(df)[7] <- "value"  ## make reshape smart     df <- cast(df, SMOKER ~ GEO) ## reshape data     row.names(df) <- df$SMOKER ## rename rows     df <- df[2:ncol(df)] ## reset df     df <- df[-4,] ## Let's only look at people who have smoked     df <- df[c("AB","BC","ON","QC")] ## and only the biggest 4 provinces (KISS)     plot(ca(df))          summary(ca(df))      Output               Principal inertias (eigenvalues):           dim    value      %   cum%   scree plot                     1      0.002523  99.9  99.9  *************************      2      3e-06000   0.1 100.0                                 3      00000000   0.0 100.0                                        -------- -----                                       Total: 0.002526 100.0                                                Rows:         name   mass  qlt  inr    k=1  cor ctr    k=2 cor ctr       1 | Crrn |  265 1000  191 |  -43 1000 191 |    1   0  43 |     2 | Dlys |  201 1000  351 |  -66 1000 351 |   -1   0  70 |     3 | Frmr |  470 1000  432 |   48 1000 432 |   -1   0  98 |     4 | Occs |   65 1000   26 |   31  964  25 |    6  36 789 |          Columns:         name   mass  qlt  inr    k=1  cor ctr    k=2 cor ctr       1 |   AB |  116 1000  146 |  -56 1000 146 |   -1   0  34 |     2 |   BC |  142 1000  775 |  118 1000 776 |   -1   0  41 |     3 |   ON |  434 1000    7 |   -6  909   6 |    2  91 540 |     4 |   QC |  308 1000   72 |  -24  994  72 |   -2   6 385 |      Looking at `summary(ca(df))` I see that nearly 100% of the inertia is described by the row profile for both modalities (Type of smoker and Province, respectively). ![CA of Smoker Types in ON, QC, AB, and BC](http://i.stack.imgur.com/ygmU7.jpg) What (I think) should be immediate takeaways are:   1. You are more likely to be a daily smoker if you live in AB, QC, or ON   2. You are more likely to be a former smoker if you live in BC    3. You are least likely to be a daily smoker if you live in BC (this fits with Canadian wide understanding of BC's "active lifestyle" culture) What could we say about occasional smokers? What would your analysis tell us through this correspondence plot and it's associated summary? Data Source: Statistics Canada, Canadian Community Health Survey (CCHS 3.1), 2005. The CANSIM table 105-0427 was an update of CANSIM table 105-0227. More current data are in CANSIM tables 105-0501 and 105-0502.


title: Distance between empirically generated distributions (in R)

text: I'm not a statistician, but I sometimes need I play around with data. I have two data sets, lists of values in the unit interval. I've plotted them as histograms, so I have an intuitive idea of how "far apart" they are. But I want to something a little more formal. My first thought was to just sum the the differences of the values in the bins, but this isn't that satisfactory. Then I thought of taking a three bin average and sum differences over these. (Apologies if I'm mangling statistics terminology) But I was thinking I'm probably reinventing the wheel, so I came here. Similar questions seem to point to "Kolmogorov Smirnov tests" or something like that. So my question is this: is this the right method to calculate how far these data sets are apart? And is there an easy way to do this in R? Ideally just `KStest(data1,data2)` or something. **Edit** To emphasise, I'm particularly interested in ways to measure how far the data are apart _directly_ rather than fitting a distribution to each and then measuring the distance between distributions. [Does that even make sense? I guess numerical calculations in R will be done by sampling from a distribution anyway...]


title: Statistics based on fractal mathematics

text: I am looking for books / textbooks on statistics based on fractal mathematics. I know it is not a very well known area and it is rather difficult to find good literature. Any suggestions are welcome (books, textbooks, online materials).


title: How to explain borderline p-values to non-stats people

text: I ran an OLS regression, and the p-value of one variable is 0.065, which is more than the threshold to consider statistically significant. However, it's not too bad either, as my data has only 300+ data points. Therefore I think the explanation needs to be taken with a bit of salt. How would you explain such a finding to non-stats people?


title: How come we can't conduct t-test and f-test in MLE?

text: Isn't the equation derived using MLE a regression? Doesn't that allow for t and f test? Please help, thanks.


title: Case weighted logistic regression

text: I'm looking at a few logistic regression issues. ("regular" and "conditional"). Ideally, I'd like to weight each of the input cases so that the glm will focus more on predicting the higher weighted cases correctly at the expense of possibly misclassifying the lower weighted cases. Surely this has been done before. Can anyone point me toward some relevant literature (Or possibly suggest a modified likelihood function.) Thanks!


title: Pivotal to estimate lambda of a exponential

text: I am studying interval estimation by the method of pivotal quantities. Let $X_1,X_2,...,X_n$ be a random sample from a p.d.f $f(x;\lambda)=\lambda e^{-\lambda x}, x>0,\lambda >0$.   I have to show that $2\lambda\sum X_i$ is a pivotal for estimating $\lambda$. Since $X $~$ exp(\lambda) $, $\sum X_i$~$Gamma (n,\lambda)$ right? so the distribution of $Q= 2\lambda\sum X_i$ has $Gamma (n,\lambda)$ distribution.This **distribution depends on $\lambda$**.Then how can I show that $2\lambda\sum X_i$ is a pivotal for estimating $\lambda$.


title: How to visualize a lot (>1k) of time spans?

text: I want to graph some data, containing of an event ID and a date (minimum required resolution is one day). every event ID has multiple dates connected to it, so the min and max date give me a time span. Another useful information would be where in time most of these events took place (like a density graph over this time span). The problem is, that i have about 1500 different events, with each more than 100 dates connected. In sum i have 1.1 million event/date pairs. (the whole data set is much larger: 200k events with an average of 10 dates) Are there any good ways to visualize this data? I tried using dot-plots, but it is more a mess and you see more a black wall than any thing else. If there is no such solution i can break down the data into smaller chunks, but for example the top 10 events (having the most dates connected) sum up to 500k tuples. If you have concrete examples in ggplot2 or D3.js it would be perfect! This is an image, how a scatter plot of start date vs end date for the top 1500 events look like. The problem is that most of the events have the same (or similar) span. ![Sample Scatter plot of the top 1500 events](http://i.stack.imgur.com/5c76Z.png) This is the same plot, but with more data. As you can see, the less common events have shorter span. ![Sample Scatter plot of the top 10000 events](http://i.stack.imgur.com/hGikC.png) Also the same data as in the last plot, but with a marker for the count of dates inside the span. As you can see there are two big blobs and the rest is quite small. ![scatter plot with count of events in span as size](http://i.stack.imgur.com/BErXb.png) now a scatter plot for all the data, you can clearly see there is a line with lots of events with only one date ![enter image description here](http://i.stack.imgur.com/UCekU.png) and as last image, this is the whole data set as density2d: ![density2d map of the whole data set](http://i.stack.imgur.com/ohBq8.png) and the same for all events with more than 10 dates: ![enter image description here](http://i.stack.imgur.com/kG2xH.png) The problem with the last two plots is, that they do not count the occurences inside the span. they just show me where a lot of time spans are. How can i put the count of events inside the time span also into this calculation?


title: Standardizing inputs for CART

text: I know I do not need to standardize the predictor variables before applying CART but would there be any adverse effects to doing it anyway? I'm comparing CART to a linear regression where I did standardize the inputs and for the sake of consistency I'd like to do the same for CART, but could it be detrimental to my results? I'm not really sure how CART goes about finding the splitting point at each node (I can't find a copy of the original paper and most of what I've read just skips over that detail) but I would think that standardizing the inputs would actually help speed up that process. Does anyone have any intuition in this regard?


title: Regression that creates $x \log(x)$ functions?

text: I just saw this website, that offers different online regression tools. For example it offers linear regression that creates a function in the form $y=ax+b$. Or it offers power regression that creates a function of the form $y=cx^a$. Or even it provides logarithmic regression that produces a function $y=c\cdot\log(x) +b$. Now, for a project I am doing, I must find functions in the form $y=c\cdot x\cdot\ln(x)+ b$. Is there any model that can do that? Also is there any tool available like the website I mentioned that offers this kind of regression?


title: A doubt about the cointegration tests

text: I'm doing Phillips-Ouliaris Cointegration test with po.test function inside tseries library. I have a simple question about the "cointegration". When i do PO test I get the p-value result but I don't get the "coefficients" of the cointegration. Is it normal? I mean, the cointegration is used to know the "relathionship" of two or more series, analyzing p-value I could understand if the series are cointegrated or not, but IF those series are cointegrated how to understand what coefficients establish their relationship? Do I have to do a linear regression AFTER the cointegration test? (to calculate the coefficients) Thank you


title: Survival Analysis—Equal follow up time

text: My question is whether my data are a survival analysis problem. I have 3 groups (2 interventions and 1 control). The 200 patients are **all** followed up for 6 months only. The follow-up is done by observing the records in a public hospital for any readmissions. I am making the assumption that a patient does not visit any other hospital. I am interested in the time (in days) to readmission of the patients after their initial discharge. The time of the initial admission was the time they entered the study. The start date is not the same for all individuals – a range of 2 years, more or less. Basically I have different start times, but the censoring time is the same for all censored patients. It is 180 days for all censored patients.   * Does this affect the survival curve or the comparison between the groups?    * Is there another analysis that I should apply that nbetter matches my data and setting? When running the Kaplan Meier curve in SPSS, I get a straight line (going down) at the 180 days ( _x_ axis). Is this OK?


title: Partial Parameter Estimation (MLE/LMMSE)

text: I have a basic question. MLE/LMSSE is introduced as follows: $$Y = H\theta + W$$ where $H$ is the linear model matrix, $W$ is measurement noise (let's assume it is normal so MLE = LMSSE). $\theta$ is the vector of parameters. The ML estimate is well-known: $(H'H)^{-1}H'\bf{Y} = \widehat{\theta}\quad\quad$ ($\bf{Y}$ is the vector of observations.) Let's say I am only interested estimating one of the parameters (say the first parameter).   Is it equal to the first entry of $\widehat{\theta}$? More generally is the ML estimate of $S\theta$ where $S$ is some low-rank matrix the same as $S\widehat{\theta}$? I do not think that is the case, but I don't know how exactly to deal with it.


title: Statistical test for measure of association not assuming monotonicity in small samples (n=6)?

text: I have two continuous variables which I have data from a physics experiment. I want to test for association between the two variables but without assuming a monotonic relationship. I also only have 6 data point each with a large error associated with it and want the test to take this into consideration. Does anyone know of a statistical test of this type?


title: Image classifier in python for few samples

text: I have 150 pictures that represent archeological signs and 5 categories to which they belong. These pictures have features like circularity, roughness and elongation that are expressed as continuous values, and other values like number of end points, number of contours, and angles between segments. I have read that SVM needs at least 100 images for every category to work well, nearest neighbour use euclidean distance metric that not consider all features in different ways, mix everythong! Alternatives? **What would be the best classifier?**


title: What does it imply when an estimate is not inside its 95% confidence interval?

text: What does it actually imply when a 95% CI does not contain an estimate (coefficient or parameter). Is there some model assumption that has not been satisfied? Or it means something else? I know when comparing two group means, if the 95% CI does not contain 0, we can conclude that two group means may be different (since m1 - m2 = 0 is not inside the 95% CI). However, I don't know how to explain why an estimated slope or intercept is not inside its 95% CI. Could it be that the estimated slope differs from the its mean too much?


title: Bicubic/bilinear interpolation in R

text: I have a data set of `x,y,z` data and I'd like to do a bicubic interpolation. x and y are spatial coordinates and z is a temperature. Below there are two images. The first one is a (gnuplot) plot of my dataset and the second one is an interpolated version (`set pm3d interpolate 10,10`). Now I'd like to do this interpolation in R but I want a matrix of values as result and not a plot. The `climates` package seemed to provide this in R, but it does not work in R 2.13 anymore. Is there another way/package to interpolate in R the way I want? ![enter image description here](http://i.stack.imgur.com/tqWCk.png) ![enter image description here](http://i.stack.imgur.com/Fyo5H.png)


title: Is it allowed to incorporate static covariates in Mixed effects models with time series data?

text: I was recently told that it is not allowed to include static covariates in mixed effect models. I can believe that this could be true, however, when I have time series data, and season is a grouping factor, then the differences between the static covariates in the different seasons should be tested, or not? If it is not allowed, how could the effect of the static covariates then be tested for? only with a repeated measures anova? sorry for this confusing question.


title: Characterizing periodicity level of signal

text: What are some good ways of quantitatively characterizing the periodicity level of an approximately periodic $f(t), t \in \mathbb{R}$ signal? I need this to tell if the output of some system is periodic or not, and to be able to decide the output of which system (out of a given set) is "more" periodic (approximates a perfectly periodic signal better). I'm interested both in purely mathematical definitions for $-\infty < t < \infty$, and in practical ways of calculating this periodicity measure for a finite and time-discretized dataset. (Note that I expect that a useful measures might give higher periodicity for longer samples of the same signal for obvious reasons---but I'd like something that works for a series of any length.) The measure should (theoretically) be the same for infinitely long perfectly periodic signals, regardless of signal shape. **EDIT:** Sample signals: http://ge.tt/91Zniy9 These signals are always between 0 and 1. I realized that the amplitude is also very relevant for my particular applications (very low amplitude signals should not be considered "periodic"), but any solution is easily modified to include this. I'm mentioning it to point out that solutions that are explicitly amplitude-dependent are also useful.


title: When the population has a known fixed size, are tables for the t statistic wrong?

text: As I understand it, Student's t will asymptotically approximate the Z statistic at sample sizes approaching infinity (see the infinite degrees of freedom row for degrees of freedom in most introductory statistics books). However, Z tests are applicable when populations are assessed in total. Thus, when the population size is finite, known, and relatively small, it seems that a t-statistic may be overly conservative. Specifically, if the population is 50 and you sampled 49, it seems that the appropriate critical value would be closer to the associated critical value for Z (given the alpha criterion) than it would be to t (given the alpha and degrees of freedom). Is it true that t is conservative in these small population cases? If so, is there a correction to the t critical values to account for cases where the population is small?


title: Counter-intuitive correlations vs. observed proportions

text: I'm confronting something pretty weird. I've got a household-level dataset, and I want to calculate the correlation between a binary variable A (taking values yes (2) or no (1)) and whether the household is urban or rural (a dummy taking values 0 or 1). I'm using Stata and the correlation I get is positive and significant at 95%. My interpretation: urban households are more likely to say yes. However, when I just look at the relative proportions of variable A for rural households and urban households, I see that urban households actually say yes _less_ often. I'm completely baffled. If you're observing a positive correlation between two dummy variables, shouldn't the observed relative proportions naturally follow that as well? I've had this happen for other variables in a completely different dataset as well. I've double- and triple-checked my Stata code, just to make sure that there's not some stupid mistake along the way, but it all seems kosher. Is this possible? What interpretation do I give it?


title: Using predict() function in R

text: After running a regression of the form `reg <\- lm(y ~ x1 + x2, data=example)` on a dataset, I can get predicted values using               predict(reg, example, interval="prediction", level=0.95)      I'm wondering what the predicted values actually refer to when I'm using the regression to predict the actual dataset. Shouldn't I obtain the original values?