title: Binary logistic regression with interaction terms

text: I have been reading several CV posts on binary logistic regression but I am still confused for my current situation. I am attempting to fit a binary logistic regression to a series of continuous and categorical variables in order to predict the mortality or the survival of animals (`qual_status`). Please see the `str` below:               > str(logit)     'data.frame':   136 obs. of  9 variables:      $ id         : Factor w/ 135 levels "01001","01002",..: 26 27 28 29 30 31 32 33 34 35 ...      $ gear       : Factor w/ 2 levels "j","sc": 2 1 1 2 1 2 1 2 2 1 ...      $ depth      : num  146 163 179 190 194 172 172 175 240 214 ...      $ length     : num  37 35 42 38 37 41 37 52 38 37 ...      $ condition  : Factor w/ 4 levels "1","2","3","4": 1 1 4 1 4 2 2 1 2 1 ...      $ in_water   : num  80 45 114 110 60 121 56 140 93 68 ...      $ in_air     : num  60 136 128 136 165 118 220 90 177 240 ...      $ delta_temp : num  8.5 8.4 8.3 8.5 8.5 8.6 8.6 8.7 8.7 8.7 ...      $ qual_status: Factor w/ 2 levels "0","1": 1 1 2 1 2 1 2 1 1 1 ...      I have no issues fitting an the following additive binary logistic regression with the `glm` function: `glm(qual_status ~ gear + depth + length + condition + in_water + in_air + delta_temp, data = logit, family = binomial)` ...but I am also interested at how these predictor variables interact with one another and possibly influence survival. However, when I attempt the following interactive binary logistic regression: `glm(qual_status ~ gear * depth * length * condition * in_water * in_air * delta_temp, data = logit, family = binomial)` I receive a warning message `"glm.fit: fitted probabilities numerically 0 or 1 occurred"`, along with missing coefficients due to singularities (NA or <2e-16 _*_ ) when I use `summary`:               Call:     glm(formula = qual_status ~ gear * depth * length * condition *          in_water * in_air * delta_temp, family = binomial, data = logit)          Deviance Residuals:        [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  0  0  0  0  0  0  0  0  0  0      [36]  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  0  0  0  0  0  0  0  0  0  0      [71]  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  0  0  0  0  0  0  0  0  0  0     [106]  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  0  0  0  0          Coefficients: (122 not defined because of singularities)                                                                 Estimate Std. Error    z value Pr(>|z|)         (Intercept)                                                1.419e+30  5.400e+22   26274077   <2e-16 ***     gearsc                                                    -1.419e+30  5.400e+22  -26274077   <2e-16 ***     depth                                                      1.396e+28  4.040e+20   34539471   <2e-16 ***     length                                                     6.807e+28  1.836e+21   37079584   <2e-16 ***     condition2                                                -3.229e+30  8.559e+22  -37727993   <2e-16 ***     condition3                                                 1.747e+31  4.636e+23   37671986   <2e-16 ***     condition4                                                 9.007e+31  2.388e+24   37724167   <2e-16 ***     in_water                                                  -4.540e+28  1.263e+21  -35935748   <2e-16 ***     in_air                                                    -4.429e+28  1.182e+21  -37470809   <2e-16 ***     delta_temp                                                -1.778e+28  3.237e+21   -5492850   <2e-16 ***     gearsc:depth                                              -1.396e+28  4.040e+20  -34539471   <2e-16 ***     gearsc:length                                             -6.807e+28  1.836e+21  -37079584   <2e-16 ***     depth:length                                              -9.293e+26  2.450e+19  -37930778   <2e-16 ***     gearsc:condition2                                          1.348e+30  3.567e+22   37809001   <2e-16 ***     gearsc:condition3                                          2.816e+30  7.495e+22   37575317   <2e-16 ***     gearsc:condition4                                                 NA         NA         NA       NA          Fitting only the continuous variables to a binary logistic regression doesn't yield any warnings or singularities but the addition of the ordinal predictor variables causes issues. Along with avoiding these warnings, is there a function/package that can handle dummy variables (I believe that is what I am looking for) in logistic regressions in `R`?


title: Log likelihood of a Markov network

text: ![enter image description here](http://i.stack.imgur.com/xCVwv.png) This is the plot of the log likelihood function from Daphne Koller's course. Two axes are $\theta_{a_1, b_1}$ and $\theta_{b_0, c_1}$. Can someone please explain for me the shape of this plot? From the equation it should be a plane isn't it?


title: How to best plot a heavy vs light tailed distribution in R

text: I want to create a plot that overlays a heavy vs light tailed distribution as an example and am trying to figure out the best way to do this. I can plot gamma distribution which is light tailed and a pareto, which is heavy but they are inherently different. As such, they are somewhat hard to compare? Any suggestions of what two distributions make a good picture? Thanks!


title: NbClust package r error

text: I try to do:               NbClust(city, diss =d, distance = NULL, min.nc = 2,      max.nc =5,method = "kmeans", index = "all", alphaBeale = 0.1)      with city is a list of number corresponding to a city and d is the dissimilarities matrix including the distance between each pair of city. I have this error:               Error in if ((!is.null(diss)) && ((distance == 1) || (distance == 2) ||  :      missing value where TRUE/FALSE needed       I try also an example on the official documentation and the program return the same error.


title: Comparing (and testing) two discrete distributions with different magnitudes

text: I am comparing authoritative survey data (large amount of observations) with data gained from a social network (very small amount of observations). Particularly, I want to compare population per district as surveyed with population per district as found in a location based social network. Example dataset:                  type variable   value     1     1      vgi    1064     2     2      vgi     873     3     3      vgi       8     4     4      vgi     246     9     1      pop 2248360     10    2      pop 3544721     11    3      pop   70934     12    4      pop 2090647      `type` is the district (1-4), `variable=vgi` denotes users found in the social network while `variable=pop` is the actual population size per distcrict. Even though the scales are completely different in magnitude, is there a way to qualitatively (e.g. with a plot) and quantitatively compare both distributions? With qualitative I mean a plot where one can easily see which district is likely _under- or overrepresented_ on social media and with quantitative I mean something like a Chi-Square-Test in order to see whether the distributions _significantly_ differ from each other. For example, one can see from the data that district `2` is underrepresented on `vgi`, or one could also say that district `1` is overrepresented on `vgi` \-- but that is the problem - what is considered over- or underrepresented?! I don't have experience with such data, thus I am asking. I was able to plot both distributions with R, but the different scales make them hard to compare - I should probably transform one of both types but I don't know how.


title: Out of 84 chances in a week, I want about 3 to be TRUE

text: I am working on an early phase software startup to help other SaaS (software as a service) companies retain new customers. We will be sending out automatic emails designed to look like they are written in a one-off personal fashion. The contents of the email will be tailored to reflect how much time they have invested in the setup of their SaaS accounts. It's our goal to send about three emails per new customer per week. Our email app sends emails during "waking" hours (8 a.m. to 8 p.m.) and it looks at the email send history for each new customer's account once per hour. So it's looking at each customer's email history 84 times per week (7*12). Upon checking the email history, it decides whether to send an email. My approach to coding this function is just to use a probability to trigger an email, with some degree of confidence that over a week the total number of sent emails is going to be about three. This will provide a lot of consistency while at the same time still appearing to be "human" in that it is random in its timing. I've tried calculating this a few times, and my best calculation                (1-(3/84))^84       has returned a probability of .04713, but I've been running a simulation for about three weeks with 2500 customer profiles and only a handful of accounts have produced 3 emails in a week. Clearly my calculations are wrong. It's been a while since I took a stat course and I'm not entirely confident I'm going to come up with the right answer :) Your thoughts are appreciated.


title: How exactly is the sum (or mean) centering constraint for splines (also w.r.t. gam from mgcv) done?

text: The data-generating-process is: $y = \text{sin}\Big(x+I(d=0)\Big) + \text{sin}\Big(x+4*I(d=1)\Big) + I(d=0)z^2 + 3I(d=1)z^2 + \mathbb{N}\left(0,1\right)$ Let $x,z$ be a sequence from $-4$ to $4$ of length $100$ and $d$ to be the corresponding factor $d\in\\{0,1\\}$. Take all possible combinations of $x,z,d$ to calculate $y$: ![enter image description here](http://i.stack.imgur.com/0X8RN.jpg) Using the (uncentered) B-spline-Basis for $x,z$ for each level of $d$ will not be feasible by the parition-of-unity-property (rows sum to 1). Such a model will not be identifiable (even without intercept). Example: _(Setting: 5 inner knot-intervals (uniformly distributed), B-Spline of degree 2, the`spline`-function is a custom one)_               # drawing the sequence     n <- 100     x <- seq(-4,4,length.out=n)     z <- seq(-4,4,length.out=n)     d <- as.factor(0:1)     data <- CJ(x=x,z=z,d=d)     set.seed(100)          # setting up the model     data[,y := sin(x+I(d==0)) + sin(x+4*I(d==1)) + I(d==0)*z^2 + 3*I(d==1)*z^2 + rnorm(n,0,1)]          # creating the uncentered B-Spline-Basis for x and z     X <- data[,spline(x,min(x),max(x),5,2,by=d,intercept=FALSE)]     > head(X)          x.1d0 x.2d0 x.3d0 x.4d0 x.5d0 x.6d0 x.7d0 x.1d1 x.2d1 x.3d1 x.4d1 x.5d1 x.6d1 x.7d1     [1,]   0.5   0.5     0     0     0     0     0   0.0   0.0     0     0     0     0     0     [2,]   0.0   0.0     0     0     0     0     0   0.5   0.5     0     0     0     0     0     [3,]   0.5   0.5     0     0     0     0     0   0.0   0.0     0     0     0     0     0          Z <- data[,spline(z,min(z),max(z),5,2,by=d)]     head(Z)              z.1d0     z.2d0      z.3d0 z.4d0 z.5d0 z.6d0 z.7d0     z.1d1     z.2d1      z.3d1 z.4d1 z.5d1 z.6d1     [1,] 0.5000000 0.5000000 0.00000000     0     0     0     0 0.0000000 0.0000000 0.00000000     0     0     0     [2,] 0.0000000 0.0000000 0.00000000     0     0     0     0 0.5000000 0.5000000 0.00000000     0     0     0     [3,] 0.4507703 0.5479543 0.00127538     0     0     0     0 0.0000000 0.0000000 0.00000000     0     0     0               z.7d1     [1,]     0     [2,]     0     [3,]     0          # lm will drop one spline-column for each factor      lm(y ~ -1+X+Z,data=data)          Call:     lm(formula = y ~ -1 + X + Z, data = data)          Coefficients:      Xx.1d0   Xx.2d0   Xx.3d0   Xx.4d0   Xx.5d0   Xx.6d0   Xx.7d0   Xx.1d1   Xx.2d1   Xx.3d1   Xx.4d1   Xx.5d1        23.510   19.912   18.860   22.177   23.080   19.794   18.727   68.572   69.185   67.693   67.082   68.642        Xx.6d1   Xx.7d1   Zz.1d0   Zz.2d0   Zz.3d0   Zz.4d0   Zz.5d0   Zz.6d0   Zz.7d0   Zz.1d1   Zz.2d1   Zz.3d1        69.159   67.496    1.381  -11.872  -19.361  -21.835  -19.698  -11.244       NA   -1.329  -38.449  -62.254        Zz.4d1   Zz.5d1   Zz.6d1   Zz.7d1       -69.993  -61.438  -39.754       NA      To overcome this problem, Wood, Generalized Additive Models: An Introduction with R, page 163-164 proposes the sum (or mean) centering constraint: $\boldsymbol{1}^T\boldsymbol{\tilde{X}_j}\boldsymbol{\tilde{\beta}_j}=0$ This can be done by reparametrization if a matrix $\boldsymbol{Z}$ is found such that $\boldsymbol{1}^T\boldsymbol{\tilde{X}_j}\boldsymbol{Z}=0$ $\boldsymbol{Z}$-matrix can be found by the QR-decomposition of the constraint matrix $\boldsymbol{C}^T = (\boldsymbol{\boldsymbol{1}^T\boldsymbol{\tilde{X}_j}})^T = \boldsymbol{\tilde{X}_j}^T\boldsymbol{1}$. _Note that $\boldsymbol{\tilde{X}_j}^T\boldsymbol{1}$ is $\boldsymbol{1}$ by the partition of unity-property._ The centered/constrained-version of my B-Spline-Matrix is:               X <- data[,spline(x,min(x),max(x),5,2,by=d,intercept=TRUE)]     head(X)              x.1d0      x.2d0      x.3d0      x.4d0      x.5d0       x.6d0     x.1d1      x.2d1      x.3d1      x.4d1     [1,] 0.2271923 -0.3225655 -0.3225655 -0.3225655 -0.2728077 -0.05790256 0.0000000  0.0000000  0.0000000  0.0000000     [2,] 0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.00000000 0.2271923 -0.3225655 -0.3225655 -0.3225655     [3,] 0.2271923 -0.3225655 -0.3225655 -0.3225655 -0.2728077 -0.05790256 0.0000000  0.0000000  0.0000000  0.0000000                    x.5d1       x.6d1     [1,]  0.0000000  0.00000000     [2,] -0.2728077 -0.05790256     [3,]  0.0000000  0.00000000          Z <- data[,spline(z,min(z),max(z),5,2,by=d,intercept=TRUE)]     head(Z)              z.1d0      z.2d0      z.3d0      z.4d0      z.5d0       z.6d0     z.1d1      z.2d1      z.3d1      z.4d1     [1,] 0.2271923 -0.3225655 -0.3225655 -0.3225655 -0.2728077 -0.05790256 0.0000000  0.0000000  0.0000000  0.0000000     [2,] 0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.00000000 0.2271923 -0.3225655 -0.3225655 -0.3225655     [3,] 0.2875283 -0.3066501 -0.3079255 -0.3079255 -0.2604260 -0.05527458 0.0000000  0.0000000  0.0000000  0.0000000                    z.5d1       z.6d1     [1,]  0.0000000  0.00000000     [2,] -0.2728077 -0.05790256     [3,]  0.0000000  0.00000000      **My question is :** Even though the fit is very similar, why do my constrained B-Spline-columns differ from what gam provides? What did I miss?               # comparing with gam from mgcv     mod.gam <- gam(y~d+s(x,bs="ps",by=d,k=7)+s(z,bs="ps",by=d,k=7),data=data)     X.gam <- model.matrix(mod.gam)     head(X.gam)       (Intercept) d1 s(x):d0.1   s(x):d0.2  s(x):d0.3  s(x):d0.4  s(x):d0.5   s(x):d0.6 s(x):d1.1   s(x):d1.2     1           1  0 0.5465301 -0.05732768 -0.2351708 -0.2259983 -0.1201207 -0.01043987 0.0000000  0.00000000     2           1  1 0.0000000  0.00000000  0.0000000  0.0000000  0.0000000  0.00000000 0.5465301 -0.05732768     3           1  0 0.5465301 -0.05732768 -0.2351708 -0.2259983 -0.1201207 -0.01043987 0.0000000  0.00000000             s(x):d1.3  s(x):d1.4  s(x):d1.5   s(x):d1.6 s(z):d0.1    s(z):d0.2  s(z):d0.3  s(z):d0.4  s(z):d0.5     1  0.0000000  0.0000000  0.0000000  0.00000000 0.5465301 -0.057327680 -0.2351708 -0.2259983 -0.1201207     2 -0.2351708 -0.2259983 -0.1201207 -0.01043987 0.0000000  0.000000000  0.0000000  0.0000000  0.0000000     3  0.0000000  0.0000000  0.0000000  0.00000000 0.5471108 -0.031559945 -0.2302910 -0.2213227 -0.1176356              s(z):d0.6 s(z):d1.1    s(z):d1.2  s(z):d1.3  s(z):d1.4  s(z):d1.5   s(z):d1.6     1 -0.01043987 0.0000000  0.000000000  0.0000000  0.0000000  0.0000000  0.00000000     2  0.00000000 0.5465301 -0.057327680 -0.2351708 -0.2259983 -0.1201207 -0.01043987     3 -0.01022388 0.0000000  0.000000000  0.0000000  0.0000000  0.0000000  0.00000000      Thw dotted line corresponds to my fit, the straight line to the gam-version ![enter image description here](http://i.stack.imgur.com/Mjptg.jpg)


title: Equations from linear regression

text: **Notation.** Let $y$, $a$, and $b$ be $n\times 1$, $p\times 1$, and $q\times1$ real vectors. Let also $X$ and $Z$ be $n\times p$ and $n \times q$ real matrices. Suppose that there is no solution, $a$, to $y = X a$. **Question**. What are the conditions on $Z$ such that $y = Xa + Zb$ has no solution for each choice of $b$? **Context**. I came across this problem in the context of linear regression. The fact that $y=Xa$ has no solution can be interpreted as "no hyperplane can perfectly fit the data". I am analysing an extension of this problem which has lead me to the need for finding something similar for "$y = Xa + Zb$ has no solution", but in this case $b$ is not fixed and can actually take any value in ${\mathbb R}^q$.


title: Does $\Pr(\text{Type I error})$ ever not equal $\alpha$ with continuous data?

text: This question: Comparing and contrasting, p-values, significance levels and type I error cites exceptions that arise when using discrete data. Are there other exceptions with continuous data as well?


title: Binary classification machine learning

text: I have   * data set (30,000) mapping people to incomes(<=some number ,>some number)    * each instance has 15 features so as age, education.  I would like some advice/pointers as to the best machine learning classifier for my task. To be implemented in Java to train. I have three main choice decision tress, navie bayes or a perception but am not sure which would best fit my problem. Any java implementations i could be pointed in the direction of would also be great. Thanks


title: What method can be used to correlate variables

text: I know a little about statistics but not much. I'm interested in learning some new methods to what seems a common problem in engineering/science. My Background: I have an engineering/science backgroun, and am pretty good with math. I have used the major mathematics software (Matlab, Mathematica, etc.), while dabbled in Python/C. I am _not_ a statistician, but I think I can pick up the fundamentals if I apply myself. I have only taken one statistics course in college. The extent of my statistics use right now is averages/standard deviation of data, with occasional linear regression. A sample problem: To answer the question on who is likely to recycle if recycling services were available and why. Perhaps there are a handful variables of interest (gender, age, geographic location, socioeconomic status, education, etc.). Let's assume I have data on a large population (probably a problem within itself). What sort of method could be used to answer this? Would this be a regression problem? Is this model development? How can I make predictions on whether someone would recycle based on the variables of interest? Extra thanks if you can point me to resources (book, website, etc.) that can help, or software I could use. Any suggestions?


title: Statistical significance of results in conversion data on group of Adwords keywords

text: I am trying to optimize my Google Adwords campaign (online advertizing on Google search), and running into the following 'statistical significance' question. I ran a campaign with many keywords (about 400 keywords), and got a total of 291 clicks (visits) and 34 conversions (clients). Now I would like to optimize my campaign, and I am focusing on the keywords statistics. I see that a set of 10 keywords _seem to_ perform much better than the rest:               Keyword group A (10 keywords - the "best"): 47 clicks, 16 conversions          Keyword group B (390 keywords - the rest):  244 clicks, 18 conversions      If I apply these data into a A/B Testing Significance calculator, I obtain a 99% confidence that A is better than B (I used a this free A/B testing significance calculator). Can I consider that group A is better than group B at converting? or am I making a wrong assumption? I am confused because all articles I have read mention A/B tests where A and B are single instances (for example just one keyword), and not group of cases (set of keywords).


title: Characteristic function of $Y=X_1X_2$, where $X_1$ and $X_2$ are standard normal

text: Let $X_1$ and $X_2$ be real-valued independent random variables with a standard normal distribution. Let $ Y = X_1X_2$. Find the characteristic function of Y. Attempt: $\phi_{Y}(t)$ = $\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty} exp(it(x_1x_2))exp(-1/2(x_1x_2)^2)dx_1dx_2$ I know that the standard normal characteristic function is: $\phi_{X_1}(t)$ = $exp(\frac{-t^2}{2})$ and likewise for $X_2$ However, I'm not sure what the resulting product is. I realize that I can use a Jacobian, but I'd like to keep it simple and see if there is a way to derive the characteristic function with something easier.


title: How can I calculate utilities for attribute levels in conjoint analysis in R?

text: I conducted a stated preference survey in which each respondent had to choose 1 set out of 3 choice sets (A, B and C), which are characterised by 4 attributes (let's say brand, color, size and keyboard yes/no) of either 2 or 3 levels. I tried to calculate the utilities/coefficients with the estimation of multinomial logit model using the mlogit package in R. It's working well to calculate the utilities for the attributes (I followed exactly the script of Yves Croissant in his paper on the mlogit package.) But I want to calculate the coefficients not for the attributes, but for the levels. I want to know what is the utility e.g. for the color "blue". How can I calculate this? Is it possible at all to calculate it with the mlogit package or R? I am grateful for every advice!


title: Pearson correlation

text: Should the sample size n be equal when we are looking for simple correlation? I mean is it OK if variable 1 has a little more or fewer observations than variable 2? I am computing correlation between two variables...the _n_ of one is a little higher than the n of other! I am looking for correlation between 2 scales (psy tests). The n of one is a little higher than the n of another. I mean not ALL the respondents who filled up one form (scale) have filled up the other. There are some (very few though) missing.


title: How to test for a critical mass - dummy variables?

text: Part of my thesis is to test where the critical mass lies. Let's assume that we want to test the influence of a certain ethnic group (x) on the performance of the whole group (members of ethnic group x and ethnic group y). The assumption/hypothesis is that if there are one or two members of group x among the whole group, the overall performance does not increase. However, if there are three or more members of group x in the group the overall performance increases. Now I have to test for this hypothesis. > Hypothesis 1a: There is no relationship between one group memeber of x and > the >performance. Hypothesis 1b: There is no relationship between two group > memebers of x and the >performance. Hypothesis 1c: There is a positive > relationship between the critical mass of group >memebers of x (at least > three) and the performance. I was thinking about using dummy variables in order to test for it. The idea was to use 0 as the reference group and then use 3 dummy variables. (Variable1: 1=1, all others = 0; Variable2: 2=1; all others=0; Variable 3: 3 or more = 1; all others=0) Then I would run a regression and see whether the Variables are significantly different from the constant. Would you recommend this procedure? I have the feeling I miss out on some important aspects. I appreciate every reccomendation and new point of view on my problem.


title: Chi-Square Test

text: my null hyp is: Age does not matter the number of male and female visitors in Park A. I perform the chi square tests with the variable of 1\. number of male and female 2\. the age. I performed the chi-square tests question 1. when there is significant association, the adjusted residue is larger than 2 or smaller than -2, but if there is any indication of a larger residue? will there be a different meaning for an adj residue of 4 and an adj residue of 2? Like the adj residue of male in younger age group is 4 while that in middle age group is 2? question 2. What if there is a negative phi value? how do I interpret the negative correlation for these two variables? Thank you so much!


title: Using anomalies to calculate trends of seasonal data

text: I commonly see people doing trend analysis of (monthly) timeseries data which show a strong inter-annual cycle following this scheme:   1. compute climatological means ("mean January", "mean February", ..., "mean December")   2. subtract climatological means from actual data, to yield an "anomaly timeseries"   3. perform linear regression on this "anomaly timeseries" _Climatological_ in this case means _multi-year average of individual months_ , e.g., an average of the 10 Januaries from 2000 to 2009. As Nick Cox points out in his comment, _anomaly_ just means deviation from a reference level; there is no implication of anything pathological or very unusual. While user31264's answer makes sense for processes where the seasonal component is truly purely additive. However, in atmospheric science we often have processes where the amplitude of the seasonal variation depends on the base level, i.e., is somewhat multiplicative. Even in these scenarios, people often use the approach I outlined above. However, I could nowhere find a rigorous statistical explanation of why this approach is actually valid. Why is the linear regression of these anomalies a reasonable solution to the regression of the original timeseries data? Can you give me any justification why this could be reasonable to do? Those people I asked mostly say "everyone's doing it" ...


title: Clarification of log interpretation

text: Let us say we have this regression $$\ln(y) = a + B_1(age) + B_2\ln(savings) + B_3\ln(income+1)$$ When carrying out the regression we obtain: $$\ln(y) = 0.3445 + 0.5(age) + 0.4556 x_1 + 0.55566 x_2$$ How would one interpret the coefficients in each case? Of particular concern is the income coefficient. An increase of income of 1% would lead to an increase of how much in $y$? This is a hypothetical example to illustrate the problem I have.


title: Bandwidth selection in multivariate semi parametric estimation in R (np package)

text: I am estimating a semi-parametric and a full non parametric model, by using the `np` package in R. I have more than 400 observations for T big years, moreover several explanatory variables will be included. This means tons of hours of work plus crashes for my machine. Therefore, I need to speed up the estimation by fixing the bandwidth to the Silverman's rule of thumb. The package has its own Silverman bandwidth, however I am trying to implement it on my own. I tried to follow `np` FAQs list, still it does not work. I am doing the following semiparametric estimation               library(np)     ## Example data.       n <- 250                    x1 <- rnorm(n)               x2 <- rnorm(n)     z1 <- rnorm(n)                ins <- rnorm(n)                y <- 1 + x1 + x2 + z1 + rnorm(n)  # dependent variable                          # setting Silverman rule of thumb bandwidth      h1 = 1.059*sd(z1)*(n^(-0.2))      h2 = 1.059*sd(ins)*(n^(-0.2))          first <- npplreg(y ~ x1 + x2|z1 + ins, bws=c(h1, h2))      I get the following error:               Error in bws[1, ] : incorrect number of dimensions      This happens also if I use the np's prebuilt Silverman rule `bw.nrd0(z1),bw.nrd0(ins)`. Note that `bw.nrd0` differs from my `h1`, this is why I want to put in my own bw However this does NOT happen if I use a full parametric form               second <- npreg(y~z1+ins,gradients=TRUE, bws=c(h1,h2))      I was not able to get what the problem is.


title: Variance Inflation Factor to Address Spatial Grouping with Binary Dependent Variable

text: I want to obtain reliable standard errors of the estimated coefficients from a regression of y on x. The observation for each individual consisted of a value of the y variable, which is binary, and a location. The latter was assigned to one of 8 geographical areas, the x values being estimated averages (from secondary data) over each of the areas. Thus the data is grouped, with x constant within each of the groups. The size of the groups varies considerably. Because of the grouping I expect that I need to apply a variance inflation factor to the OLS standard errors. Angrist & Pischke in Mostly Harmless Econometrics p 311 give this formula: $$\frac{V(\hat\beta_1)}{V_c(\hat\beta_1)} = 1+[\frac{V(n_g)}{\overline{n}} +\overline{n}-1]\rho_x\rho_e$$ Here $V(\hat\beta_1)$ is the true variance, $V_c(\hat\beta_1)$ the conventional OLS variance, $n_g$ the size of the $g$'th group, $\overline{n}$ the mean group size, $\rho_x$ the intraclass correlation of x (1 in my case) and $\rho_e$ the intraclass correlation of the error term. However, I haven't been able to find out whether this formula is valid with a binary dependent variable. Is there a formula for a variance inflation factor to address grouping in a regression with a small number of groups and a binary dependent variable?


title: Setting up a naiv tensor product B-spline example

text: I've simulated data according to $y = \text{sin}(2\cdot(4x-2))+2\cdot\text{exp}(-(16^2)(x-0.5)^2)+\epsilon$ where $\epsilon \sim \mathbb{N}(0,0.3^2)$ By evaluating the ith B-Spline of degree $k$, i.e., $\text{B}_{i,k}(x)$ for $m$ inner equidistant knot intervals recursively, I`ll get a design-matrix of dimension $n\times (m+k)$ ($q=m+k$) with B-Splines as columns. By using ordniary least squares approach it is easy to minimize $(\boldsymbol{y-X\beta})^T(\boldsymbol{y-X\beta})$ w.r.t $\boldsymbol{\beta}$. Now I wanted to set up a similar simple example for a tensor product B-spline. Here, the model would be $y_i = \sum_{j=0}^q\sum_{l=0}^p \beta_{j,l}\text{B}_{j,k_1}(x_i)\text{B}_{l,k_2}(z_i) + \epsilon_i\hspace{3ex} (1)$. W.r.t $x$, the B-Spline-Matrix $\boldsymbol{B}_x$ is of dimension $n\times q$ and w.r.t. $z$, i.e., $\boldsymbol{B}_z$ is of dimension $n\times p$. If I take the tensor product of $\boldsymbol{B}_x$ and $\boldsymbol{B}_z$, i.e., $\boldsymbol{B}_x\otimes\boldsymbol{B}_z$ this matrix will be of dimension $\Big(n^2\times (q\cdot p)\Big)$ but the vector $\boldsymbol{y}$, which will be projected onto $\boldsymbol{B}_x\otimes\boldsymbol{B}_z$, is only of length $n$? (1) suggest that, similar to the ordinary interaction terms, I take all possible products of those two B-Splines bases. But those two are column vectors? I would need to do this calculation for every pair, i.e., $q\cdot p$ different elementwise multiplications. But here the wording "tensor" would not make sense. Looking at several papers I'm still not sure how this design-matrix constructed from the tensor -product of those B-Spline matrices will look like?


title: Determine the probability density function

text: I have the following problem. Let's say `g` is a function of `A` and `R`. I run a Monte Carlo analysis to determine the pdf of `g` considering the real pdf of `R` but with `A` having a rectangular pdf. Now I have new information about the pdf of `A`. My question is: Do I need to run my Monte Carlo again or can I apply to my old pdf of `g` the newly available pdf of `A`? EDIT 1: Both old and new pdfs of `A` share the same minimum and maximum possible values, i.e. both pdfs are bounded by the same lower and upper limits. The difference is that while in the first Monte Carlo run `A` was considered to have a rectangular pdf, now with new information `A` pdf is Gaussian. It can be assumed that `A` and `R` are independent. Let the old pdf of `A` be a uniform pdf with minimum value equal to a=0 and maximum value equal to b=0.1 and the new pdf of `A` be a truncated Gaussian with equal a and b and an average value of 0.06 and a standard deviation equal to 0.01. EDIT 2: If `g` was a function of `R`, `A1`, `A2`, ..., `An` and all these random variables were independent is it also possible to update the old pdf of `g` not needing to run Monte Carlo again? Thanks!


title: How to apply GEE and how to use it for critical mass

text: sorry I am completly new to the forum and as well to using GEE and SPSS so I hope my question is not to confusing. My dependent variable is the firm performance and I have several controlling variables. First of all I am not sure how to set up my independent variable. I would like to find out whether there is a critical mass of women influencing the performance. The hypothesis is, that the performance increases when there are at least three women. However, I need to test whether this hypothesis is true. Therefore, I am looking at 8 years of panel data. That is why I was told to use the GEE method. > Hypothesis 1a: There is no relationship between one woman director and the > firm per-formance. Hypothesis 1b: There is no relationship between two women > directors and the firm performance. Hypothesis 1c: There is a positive > relationship between the critical mass of women di-rectors (at least three > women) and the firm performance. I was talking to various professors and supervisors and they all gave me contradicting opinions on how to set up the independent variable, therefore I ask again here. I hope you can help me. The second question relates to conducting a GEE analysis in SPSS with this data. I have never used this method. I know there are several tests I have to run for an OLS (normal distribution, VIF score, etc.) but I have now idea what the important tests are in the GEE method in order to make sure all assumptions hold and the model is overall making sense. It would be amazing if some of you could help me ! Laura


title: Learning multivariate techniques to analyze the gut microbiome

text: I would like to start learning some of the basics regarding data processing and analysis of microbiome data using R. Can anyone recommend a tutorial on some of the core/basic approaches that also include publicly available microbiome data sets? (Note: I understand that this would probably be highly dependent upon the questions that I want to ask about the microbiome, but surely there are fundamentals related to this type of data analysis that I could start to learn.) Thanks!


title: How to intuitively understand the relationship between difference convergences?

text: Given a series of random variables $X_1, X_2, \cdots, X_n, \cdots$ convergence in quadratic mean ⟹ convergence in probability ⟹ convergence in probability distribution How to intuitively understand this?


title: How can I estimate the density of a zero-inflated parameter in R?

text: I have a data set with lots of zeros that looks like this:               set.seed(1)     x <- c(rlnorm(100),rep(0,50))     hist(x,probability=TRUE,breaks = 25)      I would like to draw a line for its density, but the `density()` function uses a moving window that calculates negative values of x.               lines(density(x), col = 'grey')      There is a `density(... from, to)` arguments, but these seem to only truncate the calculation, not alter the window so that the density at 0 is consistent with the data as can be seen by the following plot :               lines(density(x, from = 0), col = 'black')      (if the interpolation was changed, I would expect that the black line would have higher density at 0 than the grey line) Are there alternatives to this function that would provide a better calculation of the density at zero? ![enter image description here](http://i.stack.imgur.com/1IQTU.png)


title: Why are these proper kernels and how to deduce that they are?

text: I am struggling to understand kernels and how to determine whether they are proper or not. For these examples can anyone explain why an example is proper and why another example is not. Given K1 and K2 two proper kernels. Determine which of the following formulae define proper kernels K4(x1,x2) = -K1(x1,x2)^2 K2(x1,x2) K3(x1,x2) = -K2(x1,x2) + 3k1(x1,x2) K5(x1,x2) = 15K1(x1,x2) I am particularly interested in how you get the answer, rather than stating which are proper and which are not. Thanks in advance


title: What mathematical method should be used to evaluate a questionnaire

text: This is a soft-question. I am currently evaluating a questionnaire based around social media and peoples thoughts around death for my undergrad but wanted to know what methods would one use when trying to justify the numbers to get a valid result. In that what should be the next step in evaluating the results to show that the provided response was good / bad when compared to a hypothesis. I've looked into the chi-square distribution but it seems to primarily focus on two data sets(AFAIK) such as male / female, but I was not sure if its the right path that I should follow or if there is another which would be better suited for this type of research. My question is essence: What algorithm(s) would one use when trying to collate the results from a social questionnaire when comparing to a hypothesis.


title: Testing symmetry of a distribution around its mean

text: We can test the symmetry of a distribution around $0$ by Wilcoxon sign rank test, based on its sample. But if we want to test if a distribution is symmetric around its mean, based on its sample $X_1, \dots, X_n$, is it valid to first normalize $X_i$ by the sample mean as $Y_i := X_i - \bar{X}$, and then apply Wilcoxon sign rank test to $Y_i$'s? If not, what are some ways? Thanks and regards!


title: Can anyone suggest me articles where they have used multivariate logistic regression models and explored in detail about the role of each predicor?

text: I am doing logistic regression analysis using multiple predictors for a binary outcome.I had about 10 predictors and tried to find the best model using 'glmulti' package in R. I have got a significant model with four predictors. Now I want to explore the effect of each predictor on the outcome according to the model. I thought of keeping the continuous variables at mean values and predicting probabilities for categorical predictors. Also I am trying to compute predicted probabilities for a range of continuous variable at each level of categorical predictor and plot the graph. But I need to see published medical articles where people really explored the multivariate logistic model and explained in detail about the derived model in terms of the contribution of each predictor to the outcome. Can anyone give statistical advice and also refer few articles ?


title: Using exponential smoothing to forecast irregularly spaced data in R

text: I'd like to use exponential smoothing to forecast the following data. The data is daily based. Because of some policy reasons, every $29^\text{th}$, $30^\text{th}$ and $31^\text{th}$ of each month, the data will drop to just hundreds.   1. Do I need to take out the two/three days at the end of each month since they are not following the same pattern with other days of each month?    2. If I use `Holtwinter()` to do the triple exponential smoothing, R will show error message: `time series has no or less than 2 periods.`   3. If I remove the two days at the end of each month, the time range would be `3/12/1998` to `3/28/1998`, `4/1/1998` to `4/28/1998`, `5/1/1998` to `5/28/1998`. How can I run R since the data are not in the complete month period?               18000   3/12/1998     61000   3/13/1998     59000   3/14/1998     59000   3/15/1998     66000   3/16/1998     38000   3/17/1998     37000   3/18/1998     20000   3/19/1998     72000   3/20/1998     44000   3/21/1998     37000   3/22/1998     33000   3/23/1998     28000   3/24/1998     54000   3/25/1998     24000   3/26/1998     66000   3/27/1998     52000   3/28/1998     280     3/29/1998     200     3/30/1998     400     3/31/1998     186000  4/1/1998     31000   4/2/1998     82000   4/3/1998     39000   4/4/1998     58000   4/5/1998     26000   4/6/1998     41000   4/7/1998     37000   4/8/1998     19000   4/9/1998     65000   4/10/1998     54000   4/11/1998     55000   4/12/1998     56000   4/13/1998     40000   4/14/1998     34000   4/15/1998     27000   4/16/1998     72000   4/17/1998     56000   4/18/1998     56000   4/19/1998     60000   4/20/1998     39000   4/21/1998     43000   4/22/1998     22000   4/23/1998     63000   4/24/1998     35000   4/25/1998     36000   4/26/1998     34000   4/27/1998     43000   4/28/1998     300     4/29/1998     250     4/30/1998     133000  5/1/1998     28000   5/2/1998     63000   5/3/1998     65000   5/4/1998     33000   5/5/1998     29000   5/6/1998     21000   5/7/1998     75000   5/8/1998     34000   5/9/1998     77000   5/10/1998     54000   5/11/1998     32000   5/12/1998     26000   5/13/1998     19000   5/14/1998     64000   5/15/1998     54000   5/16/1998     64000   5/17/1998     58000   5/18/1998     29000   5/19/1998     29000   5/20/1998     16000   5/21/1998     62000   5/22/1998     32000   5/23/1998     38000   5/24/1998     29000   5/25/1998     38000   5/26/1998     36000   5/27/1998     34000   5/28/1998     160     5/29/1998     150     5/30/1998


title: Weibull distribution from given mean

text: I want to pull random numbers from a Weibull distribution. The formulas for this use scale and shape parameters (both of which I have), but not a mean parameter. How can I 'build' a Weibull distribution that has a mean of my choosing? Thank you


title: Why use Chi-squared distribution to estimate the variance of normal distribution?

text: It is well known that one can use $\chi^2$ distribution to estimate the variance of normal distribution (i.e., $N(0,\sigma^2$)). Is there a particular reason of using $\chi^2$ distribution? Or can we say that it is the "most accurate" way to estimate the variance of a normal distribution? Or does this "most accurate" thing exist? If it exists and it is not $\chi^2$ distribution, then what it is? Thanks!


title: Compare two regressions

text: I try to determine India’s fertility decline between 1991 and 2001 in a multivariate regression (OLS). I have used “total fertility rate” as dependent variable and estimate the effects from six explanatory variables (see example below). I have two questions: First, How should I compare the result between 1991 and 2001? Which is the best method? Shall I investigate “interaction effects”? Or maybe “fixed- effects” and “random-effects” (seems to be a little bit complicated)? Or is it enough to simply interpret the result for each year in each column and then compare differences and similarities? Second, I would like to estimate regional effects. How shall I do this? Shall use dummy variables (like South dummy, North dummy, East dummy, West dummy). Or is a better method to use separate analysis? Which means four separate multivariate regressions for each year. Or shall I combine both methods? Example of Multivariate regression Dependent variable: Total Fertilty Rate (standard error in parenthesis)                                         INDIA 1991          INDIA 2001          Poverty                      0,014***          0,015***                                 (0,003)           (0,002)          Female literacy             -0,015***          0,0021***                                 (0,002)           (0,007)


title: How to handle data with double the number of cases than controls?

text: I have genetic SNP data with 50 controls and 120 cases. I am said that comparing few controls with more number of cases will cause spurious associations. The analysis that I want to do is just a logistic regression, where the outcome is disease (case/control) and the predictors are the genotypes (AA,AB,BB). I had a thought of comparing the 50 controls with randomly chosen 50 cases and then replicating the finding by comparing with remaining cases. Can anyone comment on this scenario.


title: The idea of making the data have a zero-mean

text: I often see people making a dimension/feature of a dataset be of a zero-mean by removing the mean from all the elements. But I never I understood why to do so? What is the effect of doing that as a preprocessing step? Does it improve the classification performance? Does it help to answer something about the dataset? Does it help when doing a visualization to understand the data?


title: PCA objective function

text: I have a quick question that I am having trouble wrapping my head around. The PCA algorithm can be formulated in terms of the correlation matrix (assume the data $X$ has already been normalized and we are only considering projection onto the first PCA). The object function can be written as: $$ \max_w (Xw)^T(Xw)\; \: \text{s.t.}$$ $$ w^Tw = 1 $$. This is fine, and we use Lagrangian multipliers to solve it. But this can be rewritten as: $$ \max_w (Xw)^T(Xw) - \lambda w^Tw $$ which is equivilant to $$ \max_w \frac{ (Xw)^T(Xw) }{w^Tw} $$ $$= \max_w \sum_{i=1}^n \text{ (distance from point $x_i$ to line $w$ ) }^2 $$ see here But this is saying to maximize the distance between point and line, and from what I've read here, this is incorrect ( it should be $\min$, not $\max$. Where is my error? Or, can someone **show me the link between maximizing variance in projected space and minimizing distance between point and line**?


title: Quantify influence in ANOVA

text: I have done the following measurements:   * Independent Variables IV1 and IV2 with a ordinal scale   * Dependent Variable 'DV' with an ordinal scale   * Tested 8 different values for IV1 and 6 different values for IV2 in a full factorial design with 5 repetition each and recorded the dependent variable. So in total 240 measurements were done. Now an ANOVA test on these measurements show that with a very high probability (> 1 - 2e-16) both independent variables contribute to the variance of 'DV'. Thats what I expected, so everything OK until here.               > aov(DV~IV1+IV2, data=d)     Terms:                     IV1                  IV2 Residuals     Sum of Squares  1721.6289       918.6018  198.5668     Deg. of Freedom         7              5       227          Residual standard error: 0.9352772      Estimated effects may be unbalanced          > summary(aov(DV~IV1+IV2, data=d))                     Df Sum Sq Mean Sq F value Pr(>F)         IV1              7 1721.6  245.95   281.2 <2e-16 ***     IV2              5  918.6  183.72   210.0 <2e-16 ***     Residuals      227  198.6    0.87                        ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1       I now want to quantify how much influence a change in IV1 and a change in IV2 has on the dependent variable DV. How can I do this? I'm not sure whether I'm asking the right question at all. A solution which came to my mind was fitting a linear model to `DV~IV1+IV2` (R syntax) an taking the coefficients as the influence. I'm not sure if this is the right way of looking at the data because I can't be sure that there is a linear dependency between DV and IV1, IV2.               > lm(DV~IV1+IV2, data=d)     Coefficients:        (Intercept)       IV1            IV2            6.2580969       0.0002853      -0.0703225


title: How to combine several time series into a useful average time series?

text: Let's assume we have four time series a, b, c and d with 10 measurments.               a(1), ..., a(10)     b(1), ..., b(10)     c(1), ..., c(10)          d(1), ..., d(10)      a, b and c are assumed to show the same trend and periodicity. The question is how can I compare d to a combination of a, b and c in order to test whether d is differing from the assumed trend and periodicity. The problem is that a, b and c have different ranges, so an average               X(i) := ( a(i) + b(i) + c(i) ) / 3      is not useful. **My question is what would be a good way to reach a meaningful combination?** Would it make sense to normalize the mean of all series a, b, c, d to 1 and then compare d to the average of a, b and c? Or would I also have to normalize the standard deviation of all four series to 1 first?


title: Is it possible to take into account a variable that isn't really specified?

text: When you run a statistic test, is it possible to take into account a variable that isn't really specified? I'll give an example. We have to see if a certain word in a sentence is collective or distributive. The words are divided into two categories: objects and subjects. We also have the two variables collective and distributive. To see whether the word in the sentence is collective or distributive, we have to look at the context of the sentence. When the context doesn't show whether it's is collective or distributive, the sentence can't be analyzed. We DO know that it has to be either one (collective or distributive), but we just don't know which one it is, because the context doesn't tell us. In this example, we want to do a Chi-square test for independence, because we want to see if collectivity and distributivity have anything to do with the position of the word (object or subject). So what do you do with those 'unspecified' sentences? I think that you can't take those sentences into account in your data analysis, because it isn't a variable. It's like Schrödinger's cat; you know that the sentence has to be collective or distributive, but as long as you don't have a good context, you can't know which one holds true. So you only take distributive and collective as variables, and the 'unspecified' sentences can't be included. Am I right or am I missing something important here?   Or is it possible to divide it like this:   * Specified -Distributive -Collective   * Unspecified So in this case you would divide the variables into: specified and unspecified. In the specified variable group you would put the distributive and the collective sentences. And in the unspecified variable group you would put all the sentences that aren't clear. And maybe even this isn't necessary. Maybe you could just have three variables: distributive, collective and unspecified. But then again, this seems really strange, because those inspecified sentences just don't seem to be a variable. Thanks in advance!


title: How to impute an ordinal variable with MICE but prevent it from taking one value?

text: I have an ordinal variable, overall_tumor_grade, that can take on values of 1, 2, 3, or X if the measurement is indeterminable. There are some NAs that I want to impute using the mice package in R, but I know that the missing values cannot be "X" because their tumor sizes are greater than 0. I want to impute overall_tumor_grade but force mice to only choose from 1, 2, or 3. Here is sample code for you to use:               df=data.frame(age=c(24,37,58,65,70,84),         overall_tumor_grade=c(1,1,2,3,'X',NA),         tumor_size=c(1.5,2.0,4.2,5.6,0,0.1))     imp=mice(df)     na_index=which(is.na(df$overall_tumor_grade))         complete(imp)$overall_tumor_grade[na_index]  #This can never be 'X'      Thank you for your help and please let me know if you need more information. * * * ## New addition @longrob suggested I temporarily remove the patients with 'X' observations, impute, then add them back in to the full dataset. Would the imputation lose power by removing all of those observations? Along longrob's suggestion, here is the work-around I have right now using a sample dataframe that is closer to what I'm really working with (several columns with missing values):               df=data.frame(age=rnorm(mean=45,sd=10,25),     overall_tumor_grade=factor(sample(c(1,2,3,'X'),25,replace=TRUE)),     tumor_size=runif(25)*10)     df[df$overall_tumor_grade=='X','tumor_size']=0 #Patients with Grade 'X' have tumorsz=0         df[sample(1:25,3),'age']=NA   ##Setting some observations to NA         df[sample(1:25,5),'overall_tumor_grade']=NA         df[sample(1:25,1),'tumor_size']=NA         ################ Imputation         imp=mice(df,meth=c('pmm',"",'pmm'))  #Suppress imputation of overall_tumor_grade         dfimp=complete(imp)         dfimp2=dfimp[dfimp$overall_tumor_grade!='X',]  #I don't want to impute grade 'X' tumors     #so I am trying to remove those observations here and then use droplevels(), but for      #a reason I can't figure out, that statement is setting all rows to NA where grade='X'      Any suggestions on how I accomplish what @longrob suggested?


title: Inputs to k-means are often normalized per-feature. Why not fully whiten the data instead?

text: We often normalize inputs to the k-means algorithm by 1) subtracting the mean on a per-feature basis and 2) dividing by the standard deviation on a per- feature basis. Some rational behind this is discussed here: Are mean normalization and feature scaling needed for k-means clustering? But it seems strange to assume that the features aren't correlated, so my question is, why don't we fully whiten the data instead? In other words, if the data has mean $\mu$ and covariance $\Sigma$, why not preprocess each sample using $ \tilde{x} = \Sigma^{-1/2} (x - \mu) $ ? The only argument I see is that this will get computationally difficult when the dimensionality of $x$ gets very large, for example $\gg 100$, but are there other reasons? Thank you.


title: What is the right name for the variant of the Kolmogorov-Smirnov statistic that retains the sign of the difference?

text: The two-sample Kolmogorov-Smirnov statistic is normally defined as $D = \max_x |A(x) - B(x)|.$ I would like to compute a variant that retains the sign of the difference between the distribution functions: $D^\prime = A(x) - B(x) \quad\text{where}\quad x = \arg\max_x |A(x) - B(x)|.$ What is the right name of $D^\prime$? Note that $D^\prime$ is different than $D^+ = \max_x [A(x) - B(x)]$ and $D^- = \max_x [B(x) - A(x)],$ which Coberly and Lewis describe as “one-sided” KS- statistics. For context, the $D^\prime$ as defined above is used by Perlman et al. to profile drug-treated cultured cells. (See page 3–4 of the supplement.) They compute $D^\prime$ independently for each image feature of the sample, then standardize the $D^\prime$-values (using the mean and standard deviation of the same measurement on mock-treated samples). Correlated profiles (of standardized $D^\prime$-values) indicate drugs that have similar biological effects.


title: mediation analysis when mediator is categorical (SPSS)

text: I want to do a mediation analysis, with the following variables:               X: Independent variable: Categorical (2 levels)     M: Mediator: Categorical (5 levels)     Y: Dependent variable: Continuous      my model: ![enter image description here](http://i.stack.imgur.com/O8QQC.png) Following MacKinnon, Lockwood, Hoffman, West, & Sheets (2002), I need to perform 2 analysis:   1. Regression analysis with X as independent and Y as dependent variable.    2. Regression analysis with X and Y as independent variables and M as dependent variable.  A full mediation would mean:   * The Regression weight X -> Y in the first analysis is significant   * The Regression weight Y -> M in the second analysis is significant    * The Regression weight X -> M in the second analysis is not significant I have 3 questions:   1. Is the above a (so not the) right thing to perform mediation analysis?    2. Is my second analysis a multinomial logistic regression?    3. Is my first analysis a GLM univariate?  For my analysis, I can only use SPSS 19 (without SEM). Reference: MacKinnon, D.P., Lockwood, C.M., Hoffman, J.M., West, S,.G., & Sheets, V. (2002). A comparison of methods to test mediation and other intervening variable effects. Psychological Methods. 3, 309 - 327.


title: bayesian logistic regression - gaussian distribution on parameters?

text: I'm trying to read this article about Bayesian logistic regression. In general, to classify instances, they use: $p(y=+1 |\beta) = \sigma(\beta^TX) $ (where $\sigma$ is obviously the sigmoid function). Now they need to find the values of $\beta$. They then say that to avoid overfitting, they'll use a prior distribution on $\beta$ specifying that each value of $\beta$ is likely to be near 0. They then impose a univariate Gaussian prior with mean =0 and variance >0 on each parameter $\beta$. (Section 3). What I don't understand is:   1. How does specifying $\beta$ to be near 0 help avoid overfitting?   2. How could one assume a Gaussian distribution on a set of _parameters_ that they're supposed to _estimate_? Any kind of help would be greatly appreciated...


title: How do I deal with asynchronous data in financial time series?

text: I have tick by tick data of two financial time series. I am trying to do online regression between the given two time series. But I am stuck due to asynchronic nature of given financial time series which causes my matrices in online regression to be singular. To deal with this, I am updating with the previous value. But this does not seem to help me much. So, how can I deal with the asynchronic nature of the given financial time series ? Any suggestion or any mathematical trick would be much appreciated.


title: how does sd of empirical log survival function and empirical log survival function

text: a sample of size 100 from an exponential distribution with parameter lambda 1. In this case, I believe that standard deviation of the empirical log survival function is `sqrt((1-e^-t)/(100*e^-t)` and log survival function is -t. when we make plots of sd of empirical log survival function and log survival function, how these two plot related?![enter image description here](http://i.stack.imgur.com/w7jYr.png) This is a plot of sd of the empirical log survival function.(x-axis:t, y-axis:sd) ![enter image description here](http://i.stack.imgur.com/dvUCz.png) This is a plot of log survival function(x-axis:t, y-axis:log(survival rate)


title: Use Random Forest model to make predictions from sensor data

text: say I have a sensor that measures temperature, pressure ++, and want to use this data to predict some quantity "A". If I use multivariate regression, I can simply implement a model of the form A=a0+a1x1+a2x2+..., and whenever I have new measurements I can use the model to make predictions. If I on the other hand make a predictive model using random forests, I'm not really sure how to use it. I've used the caret package to split my data into training and test sets, and do automatic feature selection using random forest and cross-validation. I get good predictions on the test set, but have no idea how to implement these trees to use in say a digital signal processor. In R I just use the predict() function, but this is obviously not available outside of R. This is probably a stupid questing, but it's the best I can do. Any suggestions are welcome.


title: Cophenetic distance matrix to a dendrogram

text: In hierarchical clustering procedure, a distance matrix is used to construct a dendrogram with an appropriate method of clustering. In the process of constructing a dendrogram, a **cophenetic matrix** is computed. I understand that such a cophenetic matrix is used to assess clustering consistency. How can I use a cophenetic matrix as such to plot the dendrogram ?


title: Time series with 3 factors from different sources

text: I've been asked to try to do a preliminary 'meta-analysis' of some data, but while comfortable with maths I don't have much experience with stats. I'm hoping to get some feedback with what I have attempted so far, to see if I'm on the right track. Using SPSS for the software. **Design** One dependant variable 'size' taken from different studies performed over a range of years. Most studies did not go for more than one year. Other factors:   * 'geography type' - all different types were sampled, some more than others.   * 'geography zone' - all of 4 different zones were sampled   * 'location' - many different locations were sampled, some more than others, some repeated over different years but for different studies, some only once!   * 'technique' - different sampling techniques used. **Questions**   * Is size dependant on year?   * Does geography type, geography zone or technique affect size?   * How should data be combined? Should multiple measurements from the one location in one year be averaged? **Model** I think the model should be   * Dependant variable: size    * Fixed-effects: geography type, geography zone   * Random-effects: location, technique **Tests performed so far** Linear Regression:   1. Linear regression of size vs year showed significant trend. (0.002)   2. Linear regression of size vs year, when multiple measurements from the one location for one year were averaged (one data point per location), showed not-quite-significant trent (0.060)   3. Linear regression of size vs year, when averaged location values were averaged for the entire year (one data point per year), showed no significant trend Question - which is more appropriate? ANOVA:   1. Geography type vs size: Levene statistic reported significance of .000, so ANOVA with Games-Howell option checked was performed. Between groups sig of .000   2. Geography zone vs size: Levene statistic reported significance of .007, so ANOVA with Games-Howell option checked was performed. Between groups sig of .001   3. Technique vs size: Levene statistic reported significance of .000, so ANOVA with Games-Howell option checked was performed. Between groups sig of .001 Now I think some of these might be varying together.. How to check? So what next? EDIT: MADE A MISTAKE IN THIS SECTION, have now updated Univariate GLM:   1. Fixed effects: geography type, geography zone; geo type sig (.014), zone not sig (.885), type * zone not sig (.239)   2. Fixed effects: year, geography type, geography zone; geo type sig (.002), zone not sig (.639), year sig (.000)... type * year also sig (.000)   3. Fixed effects: year, geography type, geography zone; random effects, technique; SPSS could not calculate error term.   4. Fixed effects: year, geography type, geography zone, technique; type, year and technique all significant. type * year also significant. Is this the correct approach? Is it possible that there is not enough data? When comparing means just looking at zone, there seems to be a difference. With the GLM it seems zone doesn't matter? How is this possible?


I have found a lot on the internet regarding the interpretation of random- and fixed-effects. However I could not get a source pinning down the following: **What is the mathematical difference between random- and fixed-effects?** By that I mean the mathematical formulation of the model and the way parameters are estimated.


title: Beta Regression with Zelig

text: Is it possible to use `zelig` for beta regression (as in the `betareg` package in R)? I could not find any documentation so I'm guessing it is not implemented. Is there anything that precludes its use for beta reg?


title: Dependent vs. independent samples

text: I have been searching on the net for the term dependent and independent samples, but I couldn't find proper definition, nor to conclude what exactly it is. It would be nice if anyone of you could explain it and bestow us with the knowledge.


title: How can all-positive coefficients be chosen in 2nd order linear regression to minimize the number of coefficients?

text: (this question addresses an expanded case of How can factor-levels be automatically chosen in R to maximize the number of positive coefficients in a regression model?) I am performing linear regression (using R) on data having both categorical (factor) and numeric variables, and fitting to a model having the form `y ~ (.)^2` (i.e. including all first order and second-order interaction terms). The question is: is there is a **programmatic** way of determining a coefficient vector among the set of optimal vectors $\Theta_{opt}$, where the length of the vector is minimized under the constraint that all elements of the vector are positive. There may be cases where it is impossible to find an all-positive coefficient vector in $\Theta_{opt}$, but let's assume that the particular data which is being analyzed allows for such vectors to exist. Perhaps one could start out with an initial (least-squares) optimal coefficient vector, and then manipulate this vector based on certain rules that depend on the nature of the terms which the coefficients are associated with. I can deduce some general rules for doing these manipulations, but don't know how to algorithmically perform the manipulations in a manner such that I arrive at a minimal-length parameter vector (parameters==0 don't count towards the vector's length). This may be heading the wrong direction though...?


title: Are regressions with student-t errors useless?

text: **Please see edit.** When you have data with heavy tails, doing a regression with student-t errors seems like an intuitive thing to do. While exploring this possibility, I ran into this paper: Breusch, T. S., Robertson, J. C., & Welsh, A. H. (November 01, 1997). The emperor's new clothes: a critique of the multivariate t regression model. Statistica Neerlandica, 51, 3.) (link, pdf) Which argues that the scale parameter and the degrees of freedom parameter are not identifiable with respect to each other in some sense and that because of this doing a regression with t errors doesn't do anything beyond what a standard linear regression does. > Zellner (1976) proposed a regression model in which the data vector (or the > error vector) is represented as a realization from the multivariate Student > t distribution. This model has attracted considerable attention because it > seems to broaden the usual Gaussian assumption to allow for heavier-tailed > error distributions. A number of results in the literature indicate that the > standard inference procedures for the Gaussian model remain appropriate > under the broader distributional assumption, leading to claims of robustness > of the standard methods. We show that, although mathematically the two > models are different, for purposes of statistical inference they are > indistinguishable. The empirical implications of the multivariate t model > are precisely the same as those of the Gaussian model. Hence the suggestion > of a broader distributional representation of the data is spurious, and the > claims of robustness are misleading. These conclusions are reached from both > frequentist and Bayesian perspectives. This surprises me. I don't have the mathematical sophistication to evaluate their arguments well, so I have a couple of questions: Is it true that doing regressions with t-errors is not generally useful? If they are sometimes useful, have I missunderstood the paper or is it misleading? If they are not useful, is this a well known fact? Are there other ways to account for data with heavy tails? **Edit** : Upon closer reading, of paragraph 3 and section 4, it looks like the paper below is not talking about what I was thinking of as a student-t regression (errors are independent univariate t distribuions). The errors are instead drawn from a single distribution and are not independent. If I understand correctly, this lack of independence is precisely what explains why you cannot estimate the scale and degrees of freedom independently. I guess this paper provides a list of papers to avoid reading.


title: What is the scale of prediction from models using the poisson family?

text: If I have the model:               m = glm(x~y,family=poisson)      What is the scale on which the predictions are obtained (from, e.g. `predict(m)`)? For example, I know that if `family=binomial`, the prediction scale is log-odds.


title: Methods for comparing clustering results

text: I am doing an unsupervised clustering analysis for a genomics project. This means that I do not know when a particular clustering analysis is good or not. I am running different clustering algorithms and different 'sets of features'. What I mean with different 'sets of features' is that given a data frame, I choose different combination of columns depending on its biological importance. For instance, some variables measure things at the sequence level, while others are measuring a particular cellular process or some other feature that cannot be measured at the sequence level. I am playing around with the different outputs of these sets of features, running the algorithms with _all_ the features, or ignoring some, etc . What I want is to compare the different clusters of these different runs and see if some of my objects are being clustered similarly despite lacking some sets of features. Does this make sense? Is there any recommendation on how can I do this?


title: Modeling Rare Continuous Events

text: There is extensive literature written on the modeling of rare binary events. Exact logistic regression, rare events by Gary King, and Firth method work well for binary outcomes. But what would you suggest for rare continuous outcomes? Say, I have multiple sources of income (20,000+ sources). Taken separately, each source throughout a year generates profit only on rare occasions. Each source could have 362 days of zero profit, and 3 days of positive profit. The number of profit days slightly vary from source to source. I have collected daily profit values generated by each source into one data set. It looks like pooled cross sections. This profit is my dependent variable. Independent variables associated with it are also continuous variables. Can you provide me any hints of which framework to use? (I tried tobit model that assumes left censoring, but not sure if this is a good approach.) Can I still use Firth or rare events?


title: Calculating the information gain on the features with python

text: I'm looking for a python library that computes the information gain for the features given a training matrix. Are you aware of any?


title: Microsoft Excel formula for variance

text: According to Microsoft Excel Help: > VAR uses the following formula:   >  ![Formula](http://i.stack.imgur.com/ywIpP.gif) > > where x is the sample mean AVERAGE(number1,number2,…) and n is the sample > size. Shouldn't it be n, rather than n - 1, in the denominator?


title: Approximating $Pr[n \leq X \leq m]$ for a discrete distribution

text: What's the best way to approximate $Pr[n \leq X \leq m]$ for two given integers $m,n$ when you know the mean $\mu$, variance $\sigma^2$, skewness $\gamma_1$ and excess kurtosis $\gamma_2$ of a discrete distribution $X$, and it is clear from the (non-zero) measures of shape $\gamma_1$ and $\gamma_2$ that a normal approximation is not appropriate? Ordinarily, I would use a normal approximation with integer correction... $Pr[(n - \text{½})\leq X \leq (m + \text{½})] = Pr[\frac{(n - \text{½})-\mu}{\sigma}\leq Z \leq \frac{(m + \text{½})-\mu}{\sigma}] = \Phi(\frac{(m + \text{½})-\mu}{\sigma}) - \Phi(\frac{(n - \text{½})-\mu}{\sigma})$ ...if the skewness and excess kurtosis were (closer to) 0, but that's not the case here. I have to perform multiple approximations for different discrete distributions with different values of $\gamma_1$ and $\gamma_2$. So I'm interested in finding out if there is an established a procedure that uses $\gamma_1$ and $\gamma_2$ to select a better approximation than the normal approximation.


title: How do you test an implementation of k-means?

text: **Disclaimer: I posted this question on Stackoverflow, but I thought maybe this is better suited for this platform.** How do you test your own k-means implementation for multidimensional data sets? I was thinking of running an already existing implementation (i.e., Matlab) on the data and compare the results with my algorithm. But this would require both algorithms to work more than roughly the same, and the mapping between the two results probably is no piece of cake. Do you have a better idea?


title: R getting share of users with multiple of an element

text: I have a data frame like this               User     OS     A        Windows     A        Linux     B        MacOS     C        Linux     C        FreeBSD     C        Windows     D        Windows        What I want to do is plot two types of statistics.   1. The share of users with different number of OSs. So, it would be the fraction of the users with 1 OS, with OS etc.   2. The share of users with only one OS and more than one OS. For this I tried using `tapply(OS, User, unique)`. But don't know how to go about plotting the results. I was wondering if this is the right way and what more would I need to do to to get the plots I wanted.


title: Coordinate descent for the lasso or elastic net

text: Are there any good papers or books dealing with the use of coordinate descent for L1 (lasso) and/or elastic net regularization for linear regression problems?


title: LARS vs coordinate descent for the lasso

text: What are the pros and cons of using LARS [1] versus using coordinate descent for fitting L1-regularized linear regression? I am mainly interested in performance aspects (my problems tend to have `N` in the hundreds of thousands and `p` < 20.) However, any other insights would also be appreciated. edit: Since I've posted the question, chl has kindly pointed out a paper [2] by Friedman et al where coordinate descent is shown to be considerably faster than other methods. If that's the case, should I as a practitioner simply forget about LARS in favour of coordinate descent? > [1] Efron, Bradley; Hastie, Trevor; Johnstone, Iain and Tibshirani, Robert > (2004). "Least Angle Regression". Annals of Statistics 32 (2): pp. 407–499. > [2] Jerome H. Friedman, Trevor Hastie, Rob Tibshirani, "Regularization Paths > for Generalized Linear Models via Coordinate Descent", Journal of > Statistical Software, Vol. 33, Issue 1, Feb 2010.


title: Given a control chart that shows the mean and upper/lower control limits, how do I tell if the cause of out of control points is assignable or not?

text: ![IEEE CSDA Practice Exam - Engineering Statistics - Control Chart](http://i.stack.imgur.com/Gr3f3.png) I am given 15 points. The control limits are at +/- 3 $\sigma$. Points 1, 4, 5, 6, 7, 8, 9, 10, 11, 13, and 15 fall within the control limits. Points 2, 3, 12, and 14 are outside of the control limits, with 2 being below the lower control limit and 3, 12, and 14 being above the upper control limit. How do I know if points 2, 3, 12, and 14 are out of control caused by chance causes or caused by assignable causes?


title: Multivariate skew normal

text: In the maximum likelihood estimation of Skew Normal, how does R calculate the mean? You know the formula is \begin{equation} \mu=\frac{ \sum_{i} x_{i} W(x_{i})}{\sum_{i} W(x_{i})} \end{equation}. Where \begin{equation} W(x_{i})=\frac{\phi [ \lambda \left( \frac{x_{i}-\mu}{\sigma} \right) ]}{\Phi [\lambda \left( \frac{x_{i}-\mu}{\sigma} \right)]} \end{equation} In other words, my question is: where does R get value for $\mu$, $\sigma$ and $\lambda$ when they are yet to be computed?


title: Confidence limits for a nonlinear regression

text: I've fitted some weighted data to a nonlinear model (namely a Sèrsic light profile), using least squares. That is, I have fitted a weighted nonlinear least squares model of the form: $$I_r = I_0 \exp\left[-\left(\frac{r}{\alpha}\right)^{1/n}\right]+\epsilon_r$$ where $I_0$, $α$ and $n$ are parameters. How do I then find confidence limits for each parameter, i.e. the 1,2,3 sigmas?


title: Quantile regression forecast accuracy

text: I am doing a quantile regression in R with some data and then i want to test the accuracy of the coefficients on a another data set (hold out set). But i am not sure how to go about measuring the accuracy on the hold out set in R. Is there any ideas?


title: Using [R] to connect to a Web Service

text: I have a web service, which returns a dataset, and I would like to query it from the [R] statistical package. I'd like to know how (or if) this can be done and also how to load the first datatable, of this dataset, into memory. Thanks! EDIT: It is an XML web service


title: sufficient statistic and KL-divergence: Confusion with an equation

text: I am reading a paper, which talks about minimising KL-divergence of any arbitrary distribution over a family of exponential distribution. So, given a distribution $p$, we want to compute its approximation in the exponential family. So, we have as follows where $p(x)$ is the given distribution and $q_{\theta}(x)$ is a distribution in the exponential family.: $$ \begin{split} f(\theta) &= \textrm{KL} (p \| q_{\theta}) = \left\langle \log\left(\frac{p(x)}{q_{\theta}(x)}\right) \right\rangle_{p(x)} \\\ &= \langle \log (p(x)) \rangle_{p(x)} + \langle \log(Z(\theta))\rangle_{p(x)} - \langle \theta ^{T} \phi(x) \rangle_{p(x)} \\\ &= \langle \log (p(x)) \rangle_{p(x)} + \log(Z(\theta)) - \theta ^{T} \langle \phi(x) \rangle_{p(x)} \end{split} $$ In the last line of the steps above and the last term, we can take $\theta ^{T}$ out of the expectation operator and basically treat it as a constant, which confuses me. If I understand correctly, $\theta$ are the parameters of the approximating distribution, so mean and variance if the distribution is Gaussian, for example. So, how can this be treated as a constant considering that they also depend on $x$?


title: GARCH(1,1) with insignificant constant

text: I'm fitting a GARCH(1,1) model to some data: $Y_{t} = \sigma_{t}\epsilon_{t}$ with $\epsilon_{t} \sim t(\nu)$, $\sigma_{t}^{2} = a_{0} + a_{1}Y_{t-1}^{2} + b_{1}\sigma_{t-1}^{2}$. Estimating the parameters and standard errors I get a p-value of approx. 0.26 for $a_{0}$. Now, $\hat a_{0}$ is very close to zero. I can't seem to fit the model without the constant in R, so I'm wondering whether it would be alright to proceed working with the model. Plots of the ACF/PACF of the squared standardized residuals suggest a good fit (the other estimates are significant at the 5% level (one of them barely)). Cheers!


title: Can all the p-values from drop1 models be high even when the p-value for the full lm model is low

text: I have used the R lm() function to fit a model with one numeric and ten categorical variables. The p-value reported for the regression is very low: Residual standard error: 2.459 on 320 degrees of freedom Multiple R-squared: 0.5952, Adjusted R-squared: 0.3726 F-statistic: 2.674 on 176 and 320 DF, p-value: 1.179e-14 However, when I use drop1() to test if any of the variables can be dropped, I get very high p-values (0.5 or more) for almost all the variables, and a low p-value for one of them. Is this reasonable?


title: Are features present in just one sample relevant for SVM learning?

text: I am building classification and regression SVM (RBF) models where the features for each sample are indicator values(0, 1) for a set of features exhaustively generated from all samples. There are many cases where a given feature is only present in one sample. The number of samples is greater than the number of features and I don't feel like I really need to remove features. From a mathematical standpoint (i'm still learning the underlying theory), is there a reason to keep the 'singleton' features? When they are left in, will samples containing them ever be used as support vectors?


title: Creating and interpreting Bland-Altman plot

text: Yesterday I heard of the Bland-Altman plot for the first time. I have to compare two methods of measuring blood pressure, and I need to produce a Bland-Altman plot. I am not sure if I get everything about it right, so here's what I think I know: I have two sets of data. I calculate their mean (x value) and their difference (y value) and plot it around the axis y = mean of difference. Then, I calculate the standard deviation of the difference, and plot it as "limits of agreement". This is what I do not understand - limits of what agreement? What means 95% agreement in layman's terms? Is that supposed to tell me (provided that all points of the scatter graph are between the "limits of agreement") that the methods have 95% match?


title: How to get Sphericity in R for a nested within subject design?

text: Thank you for reading. I am trying to get sphericity values for a purely within subject design. I have been unable to use ezANOVA, or Anova(). Anova works if I add a between subject factor, but I have been unable to get sphericity for a purely within subject design. Any advice? I already read the R newsletter, fox chapter appendix, EZanova, and whatever I could find online. My original ANOVA               anova(aov(resp ~ sucrose*citral, random =~1 | subject, data = p12bl, subset = exps==1))      anova(aov(resp ~ sucrose*citral, random =~1 | subject/sucrose*citral, data = p12bl, subset = exps==1))          > str(subset(p12bl, exps==1))     'data.frame':   192 obs. of  12 variables:      $ exps     : int  1 1 1 1 1 1 1 1 1 1 ...      $ Order    : int  1 1 1 1 1 1 1 1 1 1 ...      $ threshold: Factor w/ 2 levels " Suprathreshold",..: 1 1 1 1 1 1 1 1 1 1 ...      $ SET      : Factor w/ 2 levels " A"," B": 1 1 1 1 1 1 1 1 1 1 ...      $ subject  : Factor w/ 16 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...      $ stim     : chr  "S1C1" "S1C1" "S1C1" "S1C1" ...      $ resp     : num  6.01 5.63 0 2.57 6.81 ...      $ id       : int  1 2 3 4 5 6 7 8 9 10 ...      $ X1       : Factor w/ 1 level "S": 1 1 1 1 1 1 1 1 1 1 ...      $ sucrose  : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...      $ X3       : Factor w/ 1 level "C": 1 1 1 1 1 1 1 1 1 1 ...      $ citral   : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...          subset(p12b,exps==1)        exps Order       threshold SET observ S1C1 S1C2 S1C3 S1C4 S2C1 S2C2 S2C3 S2C4 S3C1 S3C2 S3C3 S3C4 S4C1 S4C2 S4C3 S4C4     1     1     1  Suprathreshold   A      1  6.0  7.1  7.5  8.6 15.0 15.4 15.0 13.1 16.9   13 13.1 16.5   24   16   21   20     2     1     1  Suprathreshold   A      2  5.6  0.8  4.0  5.6  5.6 11.3 12.9 14.5 18.5   15 12.9 14.5   24   26   29   28     3     1     1  Suprathreshold   A      3  0.0  0.0  1.7  0.0  5.0  8.4  8.4  5.0 11.7   20 18.5 16.8   29   37   37   30     4     1     1  Suprathreshold   A      4  2.6  3.3  9.1 16.3  5.4 10.0  9.6 16.8 13.5   12 22.2 23.1   19   20   22   23     5     1     1  Suprathreshold   A      5  6.8  5.3 15.4 14.5 11.5  8.3 14.5 14.2  8.9   17 11.2 15.1   24   23   19   19     6     1     1  Suprathreshold   A      6  2.6  2.8  2.6  5.2 13.4 15.6 13.7 13.0 13.7   15 16.0 18.9   22   24   25   25     7     1     1  Suprathreshold   A      7  1.3  5.8 10.2  9.8 11.9 12.3 17.7 16.7 11.4   19 19.2 21.1   16   19   18   19     8     1     1  Suprathreshold   A      8  2.0  5.6  3.9  2.0  4.9  5.2  7.5  4.9 20.2   21  8.2  9.5   30   26   32   45     9     1     1  Suprathreshold   A      9  9.4 11.3 11.7 12.1 14.7 13.8 12.6 14.9 15.2   15 15.9 13.9   17   18   15   18     10    1     1  Suprathreshold   A     10  4.5 17.8 18.5 21.6  5.8 10.9 17.0 20.2  6.6   10 17.8 18.7   12   12   16   19     11    1     1  Suprathreshold   A     11  9.8 13.0 16.1 18.0 10.5 11.6 15.4 17.3 10.1   14 15.2 16.7   13   15   15   17     12    1     1  Suprathreshold   A     12  9.6 10.4 13.3 11.3 12.1 12.6 13.6 13.6 14.9   16 15.1 16.3   16   18   18   17      Sample output               ezANOVA( data = subset(p12bl, exps==1)  , dv= .(resp), sid = .(observ), within = .(sucrose,citral), between = NULL, collapse_within = FALSE)     Note: model has only an intercept; equivalent type-III tests substituted.     $ANOVA               Effect DFn DFd  SSn  SSd     F       p p<.05   pes     1        sucrose   3  33 4953 3263 16.70 9.0e-07     * 0.603     2         citral   3  33  410  553  8.16 3.3e-04     * 0.426     3 sucrose:citral   9  99   56  791  0.77 6.4e-01       0.066      Warning messages: 1: You have removed one or more Ss from the analysis. Refactoring "observ" for ANOVA. 2: Too few Ss for Anova(), reverting to aov(). See "Warning" section of the help on ezANOVA.


title: Data representation in a better way?

text: I ran two programs `A` and `B`, and their corresponding features are `a`,`b`,`c`, and `d`. I would like to represent their features data in a better way instead of a table as shown below to explain why A is better than B or vice-versa. First, I thought about barplot representation, but feature `a` data is very small and that makes barplot representation as an improper one. Moreover, I also thought about barplot represnetatuon of A/B (normalization), but not satisfied with the resulting barplot graph. It would be great help if you share your ideas. I would like to use R for plotting the data.               prog  a   b    c    d     A    0.8 7900  70  27     B    0.3 1920 393  43      Here, the lower `d` values is better. Therefore, I would like to use the data of `a`,`b`, and `c` to explain why one program is better than other. In this case why A is better than B.


title: Comparing salaries in different countries

text: I have a number of salary figures for roles that I want to compare across various countries. Obviously each country has their own currency and living costs. A simple way of comparing would be to use the current exchange rate. Is there a more reasonable figure of "buying power" or something similar that I could use instead? Something that may incorporate living costs, etc.?


title: Is using deciles to find correlation a statistically valid approach?

text: I have a sample of 1,449 data points that are not correlated (r-squared 0.006). When analyzing the data, I discovered that by splitting the independent variable values into positive and negative groups, there seems to be a significant difference in the average of dependent variable for each group. Splitting the points into 10 bins (deciles) using the independent variable values, there seems to be a stronger correlation between the decile number and the average dependent variable values (r-squared 0.27). I don't know much about statistics so here are a few questions:   1. Is this a valid statistical approach?    2. Is there a method to find the best number of bins?   3. What is the proper term for this approach so I can Google it?   4. What are some introductory resources to learn about this approach?   5. What are some other approaches I can use to find relationships in this data? Here is the decile data for reference: https://gist.github.com/georgeu2000/81a907dc5e3b7952bc90 EDIT: Here is an image of the data: ![Industry Momentum is the independent variable, Entry Point Quality is dependent](http://i.stack.imgur.com/gr2EL.png) Industry Momentum is the independent variable, Entry Point Quality is dependent


title: PDF for ratio of Dirichlet

text: I have a random vector $\theta$, which was generated by: `X1 ~` Dirichlet($\alpha_1$)   `X2 ~` Dirichlet($\alpha_2$)   $\theta$ = element wise `X1`/`X2` What is the PDF for $\theta$? * * * Thanks at @Matthew. The author suggests approximating a Dirichlet with independent Betas. This approximation will lose the part I am most interested in. For example: If $\alpha_1$ = $\alpha_2$ = [1, 1] then:   The pdf for $\theta$ = [20, -1] is > 0 and   The pdf for $\theta$ = [20, 1] is = 0 it is unfeasible. It can't be generated through this process. If I understood it correctly the paper would assign the same pdf value to both.


title: How to calculate standard errors in OLS without inverting the X'X matrix?

text: This may seem a trivial question but I haven't found a satisfactory answer anywhere. I need to compute standard errors in a OLS regression $y = X\beta + u$ in R **from scratch**. How can I do this without inverting $X'X$ matrix? (The reason is that speed is important.) Suppose that I already have factorized the $X'X$ matric by Cholesky and already have $\hat \beta$ and $\hat \sigma$. Thanks!


title: References to papers/books that uses a kernel to smooth a discrete distribution

text: Since a kernel, such as Gaussian, is often used to smooth out the distribution of discrete points in 1D, 2D or 3D, I believe there must be some study materials or research work that have used this, and I need this for references in my work. I would appreciate if you can point me to any book/research papers/study materials that use a kernel to smooth out some distributions. Thanks


title: How to test for overdispersion in Poisson GLMM with lmer() in R?

text: I have the following model:               > model1<-lmer(aph.remain~sMFS1+sAG1+sSHDI1+sbare+season+crop       +(1|landscape),family=poisson)      ...and this is the summary output.               > summary(model1)     Generalized linear mixed model fit by the Laplace approximation      Formula: aph.remain ~ sMFS1 + sAG1 + sSHDI1 + sbare + season + crop               +      (1 | landscape)        AIC  BIC logLik deviance      4057 4088  -2019     4039     Random effects:      Groups    Name        Variance Std.Dev.      landscape (Intercept) 0.74976  0.86588      Number of obs: 239, groups: landscape, 45          Fixed effects:                   Estimate Std. Error z value Pr(>|z|)         (Intercept)  2.6613761  0.1344630  19.793  < 2e-16      sMFS1        0.3085978  0.1788322   1.726  0.08441        sAG1         0.0003141  0.1677138   0.002  0.99851         sSHDI1       0.4641420  0.1619018   2.867  0.00415      sbare        0.4133425  0.0297325  13.902  < 2e-16      seasonlate  -0.5017022  0.0272817 -18.390  < 2e-16      cropforage   0.7897194  0.0672069  11.751  < 2e-16     cropsoy      0.7661506  0.0491494  15.588  < 2e-16                     Correlation of Fixed Effects:                (Intr) sMFS1  sAG1   sSHDI1 sbare  sesnlt crpfrg     sMFS1      -0.007                                               sAG1        0.002 -0.631                                        sSHDI1      0.000  0.593 -0.405                                 sbare      -0.118 -0.003  0.007 -0.013                          seasonlate -0.036  0.006 -0.006  0.003 -0.283                   cropforage -0.168 -0.004  0.016 -0.014  0.791 -0.231            cropsoy    -0.182 -0.028  0.030 -0.001  0.404 -0.164  0.557      It is probably overdispersed, but how exactly do I calculate this? Thanks very much.


title: How to simulate multiple dependent 2x2 tables?

text: I am doing a case-control study where we have a large number, say $N$ of 2x2 tables. For a single table we have fixed column totals, say $n_{+0}$ and $n_{+1}$ such that $n_{+0} + n_{+1} = n$, the total number of patients. Given the cell probabilities $p_{00},p_{11},p_{01},p_{10}$, we can simulate a 2x2 table by generating independently $Z_0 \sim binomial(n_{+0}, \frac{p_{10}}{p_{10}+p_{00}})$ $Z_1 \sim binomial(n_{+1}, \frac{p_{11}}{p_{11}+p_{01}})$ and setting $n_{11} = Z_1, n_{10} = Z_0$, of course, $n_{01} = n_{+1}-n_{11}$ etc. My question is, how do I simulate $N$ such tables where the tables are dependent? I am deliberately making this question a little vague as I am not sure what kind of dependence between two contingency tables seem reasonable. Of course we can make any reasonable assumption, like known row totals, to make it simpler. I appreciate any help on this.


title: Derivative of order statistics

text: \--For further background into the question, one can refer to equations 2.3 and 2.6 (page 1275 of [0])-- Define: $$g_G(x)=\mbox{med}_Y|x-Y|$$ where $X$ and $Y$ are independent stochastic variables with distribution function $G$; In the above cited paper, the authors derive a closed form expression for $g_G'(x)$ when $G=\Phi$. my questions is: How does it come about that when $G=\Phi$ and $x=q$: $$g_G'(x)=\frac{\phi(x-c^{-1})-\phi(x+c^{-1})}{\phi(x+c^{-1})+\phi(x-c^{-1})}$$ where $c^{-1}=\mbox{med}_X g_{\Phi}(X)$. That part is not clear to me. reference:   * [0]:Rousseeuw, Peter J.; Croux, Christophe (December 1993), "Alternatives to the Median Absolute Deviation", Journal of the American Statistical Association (American Statistical Association) 88 (424): 1273–1283, doi:10.2307/2291267


title: Loss Matrix Equivalent with Neural Networks and random Forest

text: I'm doing classification (0,1) on a dataset for which different types of errors should be weighted differently. IE, false positives would be weighted 10 x more than false negatives. In decision trees, specifically the rpart implementation of CART in R, I can simply pass the function a loss matrix. [0,1]   [10,0] which would result in weighting errors differently and build trees based on that. for random Forest, it SHOULD be possible to do exactly the same, since it's still based on the same algorithm. With neural networks, I would think it'd still be possible, but I would have to modify the error/loss function in such a way to account for the difference in errors. I'm hoping someone will lead me in the right direction on how to go about this. Thanks.


title: Comparing two groups: one-way ANOVA vs. regression analysis with dummy

text: I am a bit puzzled with the following issue: the outcome of one-way ANOVA test shows that the mean difference of variable y between two country samples is statistically significant. However, after pooling the two samples together and running an OLS regression with y as dependent variable, some other IV and by including a country dummy, the effect of the country dummy appears to be statistically insignificant, implying that there is no country-related effect on dependent variable y. Any explanation for this outcome?


title: Statistical test Comparison of

text: I have some ordinal data (no. of exits from opiate treatment, by type (planned unplanned exits etc.). I want to know if last years figures have a statistical probability of being different from previous years? These are independent and unpaired. The null hypothesis is that there is no significant difference The totals per year are important, but I think that it's probably more important to look at the proportions of numbers in treatment per year (which obviously varies), compared to the numbers of exits. I've asked for this data, but not received it yet. I would then have to look at comparing proportions of the totals in treatment during the relevant year, is this still ordinal? as you can see I'm very green. Which test should I use?The powers that be probably want each type of exit compared by year, as well as the range per year (5years, 5 categories).


title: How to use glmnet without the intercept?

text: I know the `glmnet()` function _cannot_ exclude the intercept by users, but does anyone know how to derive the fit without intercept from the `glmnet()`?


title: Outlier detection for heavy-tailed data

text: Applying modified z-score for outlier elimination on some data (Iglewicz and Hoaglin, 1993), I discovered that a big proportion of the data (~10%) was outside the range `abs(z)>=3.5`. Further investigation showed that the data is heavy-tailed. I assumed that the Bienaymé–Chebyshev inequality would hold for the median absolute deviation MAD, too, but obviously does not.               .sum.bool  <- function(x) c('TRUE'=sum(x),'FALSE'=sum(!x),                                 'TRUE %'=round(sum(x)/length(x)*100,1), length=length(x))          rrn <- rnorm(10000)     rrt <- rt(10000,1)          # simplified z-score for demonstration purposes     mad.outlier <- function(x)abs(x-mean(x))/mad(x) > 3     sd.outlier <- function(x)abs(x-mean(x))/sd(x) > 3          rbind(mad.n=.sum.bool(mad.outlier(rrn)),           sd.n=.sum.bool(sd.outlier(rrn)),           mad.t=.sum.bool(mad.outlier(rrt)),           sd.t=.sum.bool(sd.outlier(rrt)))      On the heavy-tailed t-distribution with 1df, 14% of the data are outside `3 MADs`.                     TRUE FALSE TRUE % length     mad.n   29  9971    0.3  10000     sd.n    29  9971    0.3  10000     mad.t 1381  8619   13.8  10000     sd.t    33  9967    0.3  10000      Can someone shed light on this property of the `MAD`/z-score in the presence of heavy-tailed distributions? What are recommendations for outlier detection for heavy-tailed data?


title: Are my data zero inflated?

text: I'm building a function/package. If the users' entry data are zero inflated then the code will log-normalise them, process them, and reverse-log-normalise and bias-correct them afterwards. If they're not zero inflated it'll just process them. Does anyone know of a test for zero inflation? Everyone knows what it IS, but the only similar question went semi-off-topic and just solved the guy's problem (How to test/prove data is zero inflated?). Ideally I'd like something like this: `ifelse(iszeroinflated(data)=true, data_to_use<-log1p(data), data_to_use<-data)` But if that doesn't exist then whatever makes logical sense to be manually put in place of "iszeroinflated". What do we think? 30% of data are zeroes? 50? More? Less? And any consideration of the shape of the probability distribution other than the zeroes? One would expect lots of non-zero low numbers as well right? And only a few high ones?


title: Computationally efficient estimation of multivariate mode

text: Short version: What's the most computationally efficient method of estimating the mode of a multidimensional data set, sampled from a continuous distribution? Long version: I've got a data set that I need to estimate the mode of. The mode does not coincide with the mean or median. A sample is shown below, this is a 2D example, but an N-D solution would be better: ![enter image description here](http://i.stack.imgur.com/CpECM.png) Currently, my method is   1. Calculate kernel density estimate on a grid equal to the desired resolution of the mode    2. Look for the greatest calculated point Obviously, this calculates the KDE at a lot of non-plausible points, which is especially bad if there are a lot of data points of high dimensions or I expect good resolution on the mode. An alternate would be to use a simulated annealing, genetic algorithm, etc to find the global peak in the KDE. The question is whether there's a smarter method of performing this calculation?


title: Probability of 4 people placed in 3 rooms

text: > A family consisting of four persons—A, B, C, and D—belongs to a medical > clinic that always has a doctor at each of stations 1, 2, and 3. During a > certain week, each member of the family visits the clinic once and is > assigned at random to a station. The experiment consists of recording the > station number for each member. Suppose that any incoming individual is > equally likely to be assigned to any of the three stations irrespective of > where other individuals have been assigned. > > What is the probability that one station has two family members at it and > the others have only one? Answer: .444 I have the answers and see multiple ways to get to the answer, but it might just be a coincidence or some relation I'm not seeing. I've tried binomial method but it wasn't consistent to get other answers. I just need a clear method for these type of problems.